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Thin  sapphire  fibers,  which  find  applications  in  high  temperature  composites 
and  optical  sensors,  can  be  grown  from  the  melt  by  the  Edge-defined  Film-fed 
Gro\\Th  (EFG)  technique.  In  this  thesis,  a  thermocapillary  model  is  developed  to 
simulate  and  illustrate  the  crystal  solidification  characteristics.   The  model,  which  is 
axisymmetric  and  appropriately  scaled,  solves  the  energy  equation  in  both  the  solid 
and  liquid  phases  separately,  and  tracks  the  movement  and  evolution  of  their 
interface  explicitly  using  a  newly  developed  method  based  on  a  combined 
Lagrangian/'Eulerian  scheme. 

The  formation  of  meniscus  profiles  is  examined  utilizing  the  concept  of 
minimization  of  free  energy.  Meniscus  behavior  influenced  by  Bond  number. 
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pressurization.  aspect  ratio,  contact  angle,  and  pulling  direction  has  been  investigated 
using  both  static  and  quasi-dynamic  models  for  the  contact  condition. 

In  the  course  of  dynamic  simulation,  the  characteristics  of  the  EFG  system 
subject  to  base  temperature  perturbation  and  pull  speed  perturbation  are  studied.  To 
both  types  of  perturbation,  the  solid/melt  interface  responds  at  the  forcing  frequency. 
Due  to  the  appearance  of  multiple  time  scales  in  the  solidification  process,  both  the 
model  prediction  and  the  experimental  observation  indicate  that  as  the  forcing 
frequency  increases,  its  effects  on  crystal  growth  are  reduced.  Furthermore,  the 
interface  behaviors  yielded  by  the  static  and  dynamic  models  for  the  contact 
conditions  at  the  trij  unction  point  are  different,  indicating  the  inadequacy  of  the 
conventional  approach.  The  capability  developed  and  insight  gained  from  this  work 
can  help  improve  the  EFG  process. 
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CHAPTER  1 
INTRODUCTION 

1.1  Background  . 

Free  and  moving  boundaries  are  encountered  in  a  variety  of  engineering 
applications.   Analysis  of  such  problems  is  often  difficult  because  the  free  and 
moving  boundaries,  whose  locations  are  unknown  a  priori  and  must  be 
determined  as  a  part  of  the  solution,  make  the  problem  nonlinear.  The 
interaction  between  phase  change  and  other  thermofluid  transport  processes  adds 
additional  complexities,  making  such  problems  intrinsically  difficult  to  solve.  For 
example,  in  crystal  growth  from  the  melt,  the  coupled  heat  transfer  and  fluid 
dynamic  equations,  conjugated  with  phase  change,  soUd  phase  heat  conduction 
and  radiative  heat  exchange  with  the  ambient  surfaces,  must  be  solved 
simultaneously  in  order  to  understand  the  characteristics  of  the  process,  the 
behavior  of  the  growing  crystal  and  the  effects  of  other  manufacturing  process 
parameters  that  determine  the  quality  of  the  crystal  (Tseng  et  al.  1992,  Brown 
1988,  Viskanta  1988). 

Mathematically,  the  free  and  moving  boundary  problems  are  challenging 
and,  except  for  a  few  simple  cases,  have  to  be  solved  numerically.  Considerable 
progress  has  been  made  in  the  solution  of  problems  involving  free  and  moving 
surfaces  using  numerical  methods.,  A  central  problem  in  the  numerical  modeling 
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of  such  processes  is  the  treatment  of  heat  and  mass  flow  conditions  that  arise  at 
the  moving  solid/liquid  interface  and  the  liquid/gas  meniscus  whose  shapes  and 
locations  are  unknown  and  vary  in  time.  Many  approximation  techniques  and 
numerical  methods  have  been  developed,  as  reviewed  by  Floryan  and  Rasmussen 
(1989),  Crank  (1984),  Ozisik  (1980),  Wilson  et  al.  (1978),  Bankoff  (1964),  Carslaw 
and  Jaeger  (1959). 

The  characteristics   of  the  crystal/melt  interface  is  dictated  by  the  Stefan 
number  since  it  determines  the  growth  rate  of  the  crystal  into  the  melt.  In  order 
to  control  the  shape/quality  of  the  crystal,  the  combined  effects  of  fluid  (melt) 
motion  and  heat  transfer  in  a  variable  domain  bounded  by  capillary  (melt/gas) 
and  Stefan-type  (solid/melt)  interfaces  need  to  be  studied.  Open  research 
problems  are  of  both  mathematical  (existence,  uniqueness  of  solutions)  and 
practical  natures.  The  main  purpose  of  this  study  is  to  numerically  analyze  free 
and  moving  boundary  problems  involving  phase  change  with  special  emphasis  on 
the  edge-defmed  film-fed  growth  (EFG)  process  for  producing  thin  sapphire 
crystals. 

1.2  Description  of  Edge-defined  Film-fed  Growth  Technique 
Sapphire  fibers  for  use  in  both  optical  sensors  and  structural  composites 
have  been  produced  using  the  edge-defined  film-fed  growth  (EFG)  process 
(Kalejs  et  al.  1983a,b,  LaBelle  1980,  Sachs  1980,  Swartz  et  al.  1975),  wherein  a 
crystallographically  oriented  fibre  is  grown  from  a  meniscus  of  molten  alumina. 
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An  inert  die  is  used  to  control  the  extent  of  the  melt  and  to  help  shape  the  melt 
and  thus  the  resulting  crystal.  The  EFG  process  can  yield  continuous  lengths  of 
sapphire  fibre  with  diameters  in  the  range  of  0.05  to  0.25  mm.  The  fibre  grower, 
schematically  shown  in  Figure  1-1  (Backman  1992),  contains  three  main 
components:  the  hot  zone,  the  puller  system,  and  the  fibre  spooling  system.  The 
hot  zone  consists  of  an  inductively  heated  refractory  metal  crucible,  with  the  edge- 
defining  dies  all  enclosed  within  a  water-cooled,  environmentally  regulated 
chamber.  The  melt  is  supplied  to  each  die  tip  from  a  capillary-fed  manifold.  The 
puller  system  consists  of  a  belt  puller  and  a  fibre  guide.  The  belt  puller  utilizes  a 
double  belt  traction  mechanism  driven  by  a  precision  stepper  motor  and  is 
monitored  with  an  optical  encoder.  After  the  fibre  exits  the  puller,  it  passes 
through  several  pulleys  and  is  wound  on  spools  under  regulated  tension.  Within 
the  hot  zone,  the  fibre  growth  process  begins  with  seeding  at  the  die  tip,  which 
establishes  a  crystallographic  orientation.   The  puller  speed  is  generally  kept 
constant;  the  meniscus  dimensions  and  fibre  diameter  are  manipulated  through 
changes  in  the  induction  coil  power  level  setting  to  adjust  die  tip  temperature. 
In  practice,  both  the  die  tip  temperature  and  the  puller  speed  may 
experience  either  intentional  or  unintentional  variations  in  time.  The  unintentional 
variations  can  result  from  the  perturbations  to  the  operating  conditions  caused  by 
fluctuations  in  the  environment,  coil  power  fluctuations,  or  motor  speed 
irregularities.   The  intentional  variations  can  be  designed  to  compensate  some  of 
the  aforementioned  perturbations  to  achieve  a  controlled  growth  process  and 
uniform  crystal  properties  and  dimensions. 
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Throughout  the  development  of  the  EFG  process,  several  problems  have 
been  encountered  that  have  potentially  severe  effects  on  the  quality  of  the  grown 
crystal.  For  example,  the  shape/quality  of  the  thin  crystals  grown  from  the  melt  is 
strongly  influenced  by  the  intricate  coupling  of  heat  transfer  and  melt  flow  in  the 
EFG  system.  The  heat  transport  environment  must  yield  a  well-controlled 
temperature  field  in  both  the  solid  and  liquid  phases  so  that  the  crystallization 
rate,  the  shape  of  the  solidification  interface,  and  thermoelastic  stresses  in  the 
crystal  can  be  controlled.   A  general  goal  of  the  design  of  heat  transfer  systems 
for  any  crystal  growth  configuration  is  to  establish  a  nearly  constant  temperature 
gradient  along  the  crystal/melt  interface  and  to  control  the  cooling  rate  of  the 
crystal.  Low  dislocation  and  defect  densities  occur  when  the  temperature 
gradients  in  the  crystal  are  low  (Kalejs  et  al.  1983a,b,  Sachs  1980  and  Swartz  et  al. 
1975).  However,  because  of  the  small  dimensions  of  the  thin  crystal  sheets,  which 
are  of  the  order  of  a  few  mm,  the  rate  of  heat  losses  to  the  surroundings  along 
the  surface  of  the  crystal/melt  could  be  large.  This  usually  results  in  large 
temperature  gradients  near  the  solidification  boundary,  which  are  undesirable.  To 
overcome  this  problem,  the  shape-defining  dies  can  be  enclosed  within  a  water- 
cooled,  environmentally  regulated  chamber,  as  illustrated  in  Figure  1-1.  Active 
control  may  be  desirable  to  produce  crystals  of  uniform  cross  section  because  the 
shape  of  the  crystal  is  constrained  only  by  capillary  action. 

A  schematic  illustration  of  the  heat  transfer  mechanisms  in  the  EFG 
process  is  shown  in  Figure  1-2.  A  film  of  melt  forms  from  the  top  of  a  die  and  is 
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solidified  by  convection  and  radiative  heat  transfer  to  the  surroundings.   It  is 
important  to  control  the  growth  process  according  to  the  heat  transfer 
characteristics  in  the  system  to  obtain  a  crystal  of  definite  diameter  and  desired 
homogeneity.   In  the  melt,  convection  driven  by  mechanisms  such  as  pulling  of  the 
seed  crystal,  buoyancy  and  surface  tension  will  affect  the  heat  transfer  process  and 
hence  the  solidified  crystals.  For  thin  crystal  growth  from  the  melt  via  the  EFG 
technique  with  a  dimension  of  a  few  mm,  the  order  of  magnitude  of  convective 
heat  transfer  is  expected  to  be  relatively  small  compared  with  that  of  the 
conductive  one.  Order  of  magnitude  analysis  of  these  heat  transfer  mechanisms 
under  real  operating  conditions  are  evaluated  in  Chapter  4.  Based  on  this 
analysis,  a  mathematical  model  is  developed  for  simulating  the  dynamic 
characteristics  of  the  crystal  growth  process.  Furthermore,  the  issue  of  an 
appropriate  scaling  procedure,  which  plays  a  critical  role  in  determining  the 
computational  efficiency,  will  be  addressed. 

1.3  Scope  of  the  Present  Work 
The  purpose  of  this  study  is  to  mathematically  analyze  the  transport 
phenomena  in  the  EFG  system  so  as  to  establish  the  relationships  between  the 
operating  conditions  and  the  crystal  shape/quality.    To  develop  a  complete 
capability  of  numerical  simulation  of  free  and  moving  boundary  problems 
involving  phase  change,  several  technical  issues  need  to  be  addressed,  including 
the  assessment  of  time  stepping,  convection  and  diffusion  schemes,  meniscus 
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behavior  and  the  contact  condition  at  the  trij unction  location,  the  solution 
algorithm  for  the  thermal  field  in  both  liquid  phase  and  solid  phase,  the  front 
tracking  technique,  and  the  adaptive  gridding  to  accommodate  time-dependent 
and  irregular  domain  configurations.   In  order  to  develop  a  comprehensive  model 
to  conduct  dynamic  simulations  of  the  EFG  process,  the  following  aspects  are 
considered  in  this  thesis.  • 

First,  the  discretization  schemes  for  the  temporal  and  spatial  terms  and 
their  interaction  in  the  time-dependent  thermofluid  computations  are  investigated. 
In  order  to  clearly  understand  the  issue,  both  linear  and  nonlinear  Burgers' 
equations  are  used  as  model  problems  to  study  the  numerical  performance  of  the 
various  schemes.  The  von  Neumann  stability  and  fast  Fourier  transform 
techniques  are  utilized  to  analyze  the  stability,  accuracy,  and  dispersive  and 
diffusive  characteristics  of  the  combined  performance  of  time  stepping  and 
convection  schemes.  The  scope  of  this  aspect  of  the  study  is  more  generic  and  is 
intended  to  provide  some  general  guidelines  for  helping  devise  suitable 
computational  strategies  for  growing  crystals  of  large  sizes,  where  convection 
effects  are  more  pronounced  and  may  vary  spatially. 

Second,  the  physics  and  appropriate  treatment  of  the  solid-liquid-gas 
trijunction  point,  which  is  not  well  understood  at  the  present  time,  is  explored. 
Consideration  of  force  equilibrium  leads  to  a  stress  singularity  at  the  contact  line. 
The  modeling  of  static  contact  lines,  where  a  contact  angle  is  specified  at  the 
trijunction  point,  is  justified  only  in  the  liquid-gas-solid  system  under  the 
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assumption  that  the  solid  wall  is  smooth  and  homogeneous.    The  extension  of  this 
concept  to  dynamic  contact  lines  is  desirable,  to  relate  the  instantaneous  contact 
angle  to  the  velocity  of  contact  line  motion.  The  attempts  to  determine  this 
function  from  physical  experiments  have  not  been  very  successful  and  this  remains 
an  open  area  of  research.   Relatively  detailed  discussion  will  be  given  later,  for 
both  the  static  and  dynamic  contact  conditions.   A  simplified  dynamic  model, 
taking  into  account  the  effect  of  the  direction  and  speed  of  contact  point 
movement,  will  also  be  proposed  and  incorporated  into  the  present  model.  Based 
on  both  the  static  and  dynamic  contact  conditions,  meniscus  characteristics  and 
solidification  behavior  in  the  EFG  process  will  be  examined  and  compared. 

The  third  aspect  of  this  study  is  to  simultaneously  solve  the  equations  of 
heat  transfer  in  both  solid  and  liquid  phases,  liquid/solid  interface  movement,  and 
the  dynamics  of  meniscus  and  contact  lines  using  a  curvilinear  coordinate  system. 
Body-fitted  coordinates  are  used  for  mapping  the  irregular  shape  of  the  timewise- 
varying  solid/liquid  and  liquid/gas  interfaces.   A  new  Lagrangian  approach  is 
developed  for  tracking  the  phase  front,  and  a  solution  procedure  is  used  for 
solving  simultaneously  the  thermal  fields  in  both  liquid  and  solid  phases.  This 
yields  a  complete  EFG  numerical  simulation. 

Finally,  thermal  sensitivity  and  stability  of  the  EFG  for  sapphire  fibre 
production  is  numerically  studied.  The  dynamic  behavior  of  the  EFG  system  with 
respect  to  an  external  disturbance,  e.g., varying  the  base  temperature  or  pull 
speed,  is  investigated.   The  proposed  dynamic  model  is  incorporated  to  study  the 
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dynamic  EFG  characteristics.    Based  on  this  model,  the  variation  of  radius  of  the 
grown  crystal  is  compared  with  the  static  wetting  condition,  which  is 
conventionally  employed  in  numerical  simulations  of  crystal  growth  (Crowley 
1983,  Ettouney  and  Brown  1983).  The  relationship  between  response  and  forcing 
frequency,  in  terms  of  magnitude  and  phase  shift,  is  identified. 

1.4  Outline 

The  present  study  therefore  can  be  divided  into  several  parts: 
discretization  schemes  for  temporal  and  spatial  terms  and  their  interaction, 
meniscus  and  contact  point  behaviors,  front  tracking  and  solution  techniques  for 
phase  change  problems,  and  dynamic  simulation  of  the  EFG  process. 

In  Chapter  2,  the  time  stepping  schemes,  the  convection  schemes,  and  their 
interaction  in  unsteady  thermofluid  computations  are  presented.    Both  linear  and 
nonlinear  Burgers'  equations  are  used  as  model  problems  to  study  the  stability, 
accuracy,  and  dispersive  and  diffusive  characteristics  of  each  scheme.  With  the 
aid  of  a  stability  analysis,  the  overall  performance  of  the  scheme  is  discussed. 

The  shape  of  the  melt/gas  meniscus  that  connects  the  growing  crystal  with 
the  die,  which  is  governed  by  the  Laplace- Young  equation,   is  critical  in 
determining  the  properties  of  the  solidified  crystal.  In  Chapter  3,  the  meniscus 
characteristics  with  application  to  the  EFG  process  are  studied.  Both  static  and 
dynamic  contact  line  motions  are  investigated.   Quasi-equilibrium  dynamics  of  the 
meniscus  is  simulated  using  a  simplified  hysteresis  model  for  the  contact  angle  at 
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the  trij  unction  point.  The  implication  of  the  variation  of  the  meniscus  shape  on 
the  crystal  growth  characteristics  is  explored  in  detail. 

In  Chapter  4,  a  general  treatment  of  the  interface  movement  is  designed, 
which  is  not  limited  to  a  single-valued  or  isothermal  surface  employed  by  other 
investigators  (Lacroix  1990,  Thompson  and  Szekely  1989,  Ettouney  and  Brown 
1983,  Crowley  1983,  Sparrow  et  al.  1977).  The  development  of  the  numerical 
technique  is  reported  in  Shyyet  al.  (1992b),  which  demonstrates  that  the  method 
is  robust  and  can  handle  both  isothermal  and  nonisothermal  phase  change 
problems.  The  technique  is  combined  with  a  moving  gridding  procedure, 
determining  the  thermal  field  in  both  the  solid  and  liquid  phases  and  the 
evolution  of  both  the  solid/melt  interface  and  the  melt/gas  meniscus  to  predict 
the  solidification  characteristics  of  sapphire  fiber  growth  via  the  EFG  technique 
for  different  Stefan  numbers.  Two  different  scaling  procedures  for 
nondimensionalizing  the  governing  equations  are  examined,  and  the 
computational  characteristics  and  relative  efficiency  of  these  two  procedures  are 
compared. 

In  Chapter  5,  the  sensitivity  and  dynamic  characteristics  of  the  EFG  system 
subject  to  perturbations  of  the  operating  environment  are  explored.  The  dynamic 
behavior  of  the  EFG  system  with  respect  to  the  external  disturbance,  i.e.,  variations 
in  either  base  temperature  at  the  die  tip  or  pull  speed,  is  investigated.   A  dynamic 
model  is  also  incorporated  into  the  system  and  compared  with  the  conventional 
static  model  for  the  contact  condition  at  the  trijunction  point.  In  Chapter  6,  a 
summary  and  some  suggestions  for  further  work  are  given. 
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Overall,  the  contributions  of  this  thesis  are  to  develop  and  utilize  advanced 
and  original  numerical  techniques  for  simulating  the  EFG  process,  to  examine  and 
assess  the  contact  conditions  relevant  to  the  present  physical  problem,  and  to 
construct  a  predictive  tool  incorporating  both  the  physical  models  and  the 
numerical  techniques.   The  developed  model  system  is  envisaged  as  a 
computational  tool,  to  make  detailed  investigations  of  the  dynamic  characteristics 
of  the  EFG  processes. 
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Figure  1-2  Schematic  illustration  of  heat  transfer  mechanisms 
in  EFG  process. 


CHAPTER  2 

INTERACTION  OF  TIME  STEPPING  AND  CONVECTION  SCHEMES 
FOR  UNSTEADY  FLOW  SIMULATION 

2.1  Introduction 

The  numerical  solutions  of  heat  transfer,  fluid  flow,  and  other  related  transport 
problems  can  be  obtained  by  solving  the  laws  governing  these  processes  in  discrete 
forms.  By  using  appropriate  numerical  techniques  to  discretize  the  spatial 
derivatives,  the  original  continuous  problems  are  reduced  to  a  system  of  ordinary 
differential  equations  in  time.  A  time  stepping  scheme  can  then  be  used  to  advance 
the  solution  in  time  from  given  initial  and  boundary  conditions.  Categorically,  two 
kinds  of  time  stepping  methods  can  be  defined,  namely,  explicit  and  implicit.  The 
explicit  schemes  provide  a  noniterative  "marching"  procedure  for  obtaining  the 
solution  of  each  nodal  point  at  present  time  in  terms  of  the  known  preceding  and 
boundary  values.  Explicit  schemes  limit  the  time  step  size  due  to  the  stability  as  well 
as  the  accuracy  requirements.  Implicit  procedures,  on  the  other  hand,  generally 
involve  the  simultaneous  solution  of  all  unknowns  in  terms  of  both  the  known 
preceding  and  the  unknown  present  values  at  other  nodes.  The  implicit  methods  are 
usually  more  stable,  and  their  restriction  in  time  step  size  is  mainly  due  to  the 
accuracy  consideration. 
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Besides  the  time  stepping  methods,  the  convection  schemes  play  an  important 
role  in  both  steady  and  unsteady  flow  computations.  Discretization  of  the  convection 
terms  is  one  of  the  principal  difficulties  in  computing  the  complex  fluid  flows.  The 
first-order  upwind  scheme  used  to  be  a  popular  choice  due  to  its  superior  numerical 
stability  as  well  as  the  wiggles-free  characteristics.  However,  the  resulting  accuracy 
yielded  by  this  scheme  is  often  found  to  be  inadequate  in  view  of  its  first-order 
truncation  error  and  strong  numerical,  as  opposed  to  physical,  diffusion  effects  (Shyy 
1985,Roache  1972). 

Many  alternative  convection  schemes  have  been  studied,  as  reviewed  in  Hirsch 
(1990),  Fletcher  (1988),  Correa  and  Shyy  (1987),  and  Warming  et  al.  (1973),  and  the 
findings  for  the  same  schemes  are  not  necessarily  consistent  among  different  studies; 
they  are  obviously  problem  as  well  as  implementation  dependent.  Furthermore,  it 
is  known  that  the  performance  of  a  given  convection  scheme  can  be  quite  different 
between  the  steady  and  unsteady  flow  problems  (Shyy  1984),  hence,  it  is  important 
to  address  the  issues  related  to  the  interactions  between  the  time  stepping  and  the 
convection  schemes  for  unsteady  flow  calculations.  In  this  chapter,  assessment  of  the 
relative  merits  of  the  various  possible  choices  for  solving  the  convection-dominated 
transport  problems  is  made  through  a  series  of  well  defined  cases. 

The  viewpoint  taken  here  is  that  a  numerical  scheme  preferably  is  of  high 
order  of  accuracy  in  truncation  error  not  only  so  that  it  can  perform  well  for 
problems  of  smooth  solution  profiles,  but,  equally  importantly,  so  that  no  over/under- 
shoots or  spurious  oscillations  that  violate  the  conservation  laws  appear,  in  order  to 
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be  able  to  honor  physical  realizabiiity.  In  the  latter  regard,  the  concept  of  total 
variation  diminishing  (TVD)  has  been  found  very  helpful  for  constructing  an 
appropriate  convection  scheme  (Hirsch  1990).  It  is  known  that,  with  the  use  of  the 
forward  Euler  scheme,  the  only  convection  scheme  that  can  guarantee  the  TVD 
property  across  a  sharp  gradient  or  a  discontinuity  is  first-order  (Roe  1986).  Here, 
through  pragmatic  testing  and  basic  analysis,  the  intent  is  to  gain  insights  of  the 
performance  of  the  schemes  considered  without  insisting  on  the  satisfaction  of  TVD. 
The  reasons  are  twofold.  No  known  technique  is  available  to  enforce  TVD  for  a 
uniformly  high  order  time  stepping  scheme;  furthermore,  for  problems  with 
substantial  source  terms,  the  concept  of  TVD  does  not  apply  (Hirsch  1990). 

Three  time  stepping  schemes,  the  first-order  backward  Euler,  the  second-order 
Crank-Nicolson.  and  the  third-order  Adams-Bashforth/Adams-Moulton  (A-B/A-M) 
predictor-corrector  schemes,  are  considered.  In  conjunction,  four  convection 
schemes,  first-order  upwind,  second-order  central  differencing,  second-order  upwind, 
and  QUICK  (Leonard  1979),  are  examined  using  both  the  linear  and  nonlinear 
Burgers  equations.  The  goal  is  to  understand  the  dispersive  and  dissipative 
characteristics  exhibited  by  each  scheme  and  to  offer  guidance  for  utilizing  these 
features  to  help  improve  the  overall  numerical  accuracy.  The  Burgers  equation,  of 
the  linear  as  well  as  nonlinear  forms,  is  used  as  the  test  problem.  Furthermore,  the 
von  Neumann  stability  analysis  and  the  FFT  spectral  analysis  will  also  be  used  to 
help  understand  the  characteristics  of  each  formulation. 
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The  results  to  be  presented  in  the  following  are  a  part  of  the  attempt  of 
developing  the  necessary  models  for  predicting  the  crystal  growth  processes  with 
different  physical  scales  and  operating  conditions.  As  will  become  clear  later,  the 
convective  effect  is  not  dominant  in  the  present  EFG  processes;  however,  the  work 
conducted  here  will  become  essential  for  any  system  of  substantial  convective  effects. 
Assessment  of  the  time  stepping  schemes  for  phase  change  problem,  which  is  critical 
for  the  problem  under  study,  is  presented  in  section  4-3. 

2.2  Model  Problem  Analysis 
The  model  problem  considered  is  the  unsteady  one-dimensional  Burgers 
equation 

^^^{u4>)=—{i>—4>),  0<x<  l,r>0,  v=comtant  >Q  (2.1) 
at     dx  dx  dx 

where  if  u  is  prescribed  constant,  the  equation  is  linear,  while  if  u  is  taken  to  be  a 
function  of  (f),  then  it  is  nonlinear.    Both  cases  are  adopted  here  since  they  are 
relevant  to  practical  applications  under  different  circumstances.  The  boundary  and 
initial  conditions  are  given  in  the  following  : 
(I)  Linear  Problem  (  u  =  1  ) 
(i)  boundary  conditions 

0(0,0  =  1 
<t>(l,t)  =  0 
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(ii)  initial  condition 

<t>(x,0)  = 

(II)  Nonlinear  Problem  (  u  =  <^  ) 
(i)  boundary  condition 


1  ,  for  0<x<0.2 
0  ,  otherwise 


0(0,0  =  0 


0(1,0  =   J  

l+(I)^exp(J-) 


(ii)  initial  condition 

<t>ix,l)  = 


l+(-)^exp(£!) 


where  to  =  exp(l/8v). 

The  corresponding  analytic  solution  (Lohar  and  Jain  1981,  Benton  and 
Platzman  1972)  is 


X 


Hx,t)  =   i   (2 

IM— )'exp(£l) 

In  what  follows,  the  second-order  central  differencing  scheme  is  always 
adopted  to  approximate  the  diffusion  term,  i.e.. 


(2.3) 


where  T,,  ,  of  the  order  of  0[(Ax)^],  is  the  truncation  error. 

As  to  the  convection  term,  the  following  schemes  are  considered  and  tested: 
(a)  first-order  upwind 


d(u4>) 
dx 


(u<t>)^  -  (M<^),., 
Ax 

Ax 


+  r  ,  foru>0 
+  r  .      foru  <0 


where  T„  of  order  of  0[(Ax)],  is  the  truncation  error, 
(b)  second-order  central  differencing 

d(u({>)  I   ^  -  iu<i>\.,  ^ 


dx 


2Ax 


where  T,  is  0[(Ax)^] 
(c)  second-order  upwind 


(2.4) 


(2.5) 


d(u<l>)  ■  _. 
dx 


-L[  3(M0), -4(M0),.,  +(M<^)...3  ]  +  r  ,  foru>0 

ZAX  (^.O) 


1 


2Aa: 


where  T,  is  0[(Ax)^] 
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(d)  QUICK 


diu4>)  I  ^. 
dx 


J_[ -{uci> ) . +7(u<^. ) . -3(M<;.) .  -3(M<^), +7  ,/or  u  <  0 
8Ax 


where     is  also  0[(Ax)^]. 

When  the  first-order  backward  Euler  method  for  the  time  derivative  is  used, 
the  discretized  form  of  Eq.  (2.1)  can  be  expressed  as  follows: 

(a)  first-order  upwind 

where  P  =  u  Ax/p  =  cell  Peclet  number, 

C  =  u  At/ Ax  =  Courant  number 
and  the  superscript  "o"  represents  the  known  value. 

(b)  second-order  central  differencing 


(c)  second-order  upwind 


20 


(d)  QUICK 


(2.11) 


The  Crank-Nicolson  scheme  is  a  second-order  accurate  method  and  usually 
is  described  as  unconditionally  stable  (Fletcher  1988,  Anderson  et  al.  1984). 
However,  as  demonstrated  later,  this  does  not  mean  that  the  solutions  are  free  of 
spurious  oscillations.  If  the  convection  and  diffusion  terms  are  replaced  by  the 
appropriate  discretization  schemes  at  both  the  current  and  previous  time  levels,  the 
corresponding  discretized  equations  result: 
(a)  first-order  upwind 


(2.12) 


(b)  second-order  central  differencing 


(Ul)  0,..  -(2-2Z)  0,.  .(1-^)  0,.^,  = 


(2.13) 


(c)  second-order  upwind 
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(2.14) 


(d)  QUICK 


(2.15) 


The  third-order  predictor-corrector  scheme,  based  on  the  Adams-Bashforth 
scheme  as  the  predictor  step  and  the  Adams-Moulton  scheme  as  the  corrector  step, 
is  explicit;  it  involves  successive  computations  at  each  node  but  does  not  need  to 
solve  the  whole  domain  simultaneously. 

If  Eq.  (2.1)  is  rewritten  in  the  following  manner, 


at         ox  dx  dx 


(2.16) 


then  the  present  predictor-corrector  algorithm  results  : 


Adams-Bashforth  method  : 


(2.17) 


Adams-Moulton  method 
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^c'l)  =       +  ^(  5/c*i)-  +  8/w  )  (2.18) 

where  (n+1)*  designates  the  value  obtained  from  the  predictor  step.  It  is 
noted  that  the  appropriate  convection  and  diffusion  schemes  can  be  used  to  represent 
the  function  f(x,<^).  The  present  A-B/A-M  scheme  needs  a  starting  procedure  since 
it  requires  information  from  the  previous  two  nodes.  Because  the  solutions  close  to 
the  left  boundary  are  smooth,  it  is  found  that  no  discernible  difference  exists  between 
using  the  Crank-Nicolson  and  the  backward  Euler  schemes  as  the  starting  procedure. 

2.3  Results  and  Discussions 
The  viewpoint  taken  here  is  that  a  numerical  scheme  preferably  is  of  high 
order  of  accuracy  in  truncation  error  not  only  so  that  it  can  perform  well  for 
problems  of  smooth  solution  profiles,  but  equally  importantly,  so  that  no  over/under- 
shoots  or  spurious  oscillations  that  violate  the  conservation  laws  appear,  in  order  to 
be  able  to  honor  physical  realizability.  In  the  latter  regard,  the  concept  of  TVD  has 
been  found  very  helpful  for  constructing,  as  well  as  assuming,  the  convection  scheme 
Hirsch  (1990).  It  is  known  that,  with  the  use  of  the  forward  Euler  scheme,  the  only 
convection  scheme  that  can  guarantee  the  TVD  property  is  first  order.  Here, 
through  pragmatic  testing  and  basic  analysis,  the  intent  is  to  gain  insights  of  the 
performance  of  the  schemes  considered  without  insisting  on  the  satisfaction  of  TVD. 
The  reasons  are  twofold.  No  known  technique  is  available  to  enforce  TVD  for  a 


23 

high  order  time  stepping  scheme;  furthermore,  for  problems  with  substantial  source 
terms,  the  TVD  scheme  does  not  work  well  anyway  (Hirsch  1990). 

The  results  of  the  linear  and  nonlinear  Burgers  equations  are  presented  in  the 
following  sections. 

2.3.1  Linear  Burgers  Equation 

The  linear  Burgers  equation  is  chosen  because  it  is  relevant  to  many  heat 
transfer  applications.  It  is  known  that  at  steady  state,  both  first-order  and  second- 
order  upwind  solutions  are  wiggles-free  while  neither  second-order  central 
differencing  nor  QUICK  is  wiggles-free  (Shyy  1985).  Their  relative  merits  in  the 
context  of  the  unsteady  problems  are  assessed  here.  The  time  stepping  schemes 
employed  here  are  forward  and  backward  Euler  schemes.  The  A-B/A-M  scheme  will 
be  included  later  in  the  nonlinear  case. 

For  the  linear  problem,  with  C  =  1  and  P  =  10,  either  21  or  81  nodal  points 
are  utilized  in  the  computations.  The  results  with  21  nodes  are  first  presented  and 
compared  in  Figures  2-1  and  2-2  for  the  two  time  stepping  schemes  at  t  =  0.3,0.8, 
and  2.  Using  the  backward  Euler  scheme,  as  shown  in  Figure  2-1,  both  the  first- 
order  and  second-order  upwind  solutions  are  wiggles-free,  but  the  first-order  upwind 
is  a  little  more  dissipative.  For  the  steady  state  solutions  of  Eq.  (2.1),  it  is  well 
known  that  the  central  differencing  scheme  exhibits  oscillations  if  the  cell  Peclet 
number,  P,  exceeds  2  (Fletcher  1988  and  Shyy  1985).  Figure  2-1  shows  that  this 
scheme  also  produces  nonphysical  oscillations  behind  the  sharp  front  for  the  unsteady 
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equation.  The  QUICK  scheme,  which  combines  the  characteristics  of  the  central 
differencing  and  upwind  schemes,  also  exhibits  oscillations  in  the  regions  of  sharp 
gradient;  its  solution  profiles  lie  between  those  of  the  central  differencing  and 
upwind  schemes. 

At  t  =  2,  as  shown  in  Figure  2-1  (c),  the  solutions  have  approached  steady 
state.  Both  the  first-  and  second-order  upwind  schemes  converge  to  virtually  the 
same  solution.  Even  though  the  first-order  upwind  scheme  works  nearly  as  well  as 
the  second-order  upwind  convection  scheme  for  the  one-dimensional  steady  state 
equation,  when  combined  with  the  backward  Euler  time  stepping  scheme,  it  is  less 
satisfactory  for  the  unsteady  equation  in  view  of  the  stronger  numerical  smearing 
effects  depicted  in  Figures  2-1  (a)  and  (b).  For  the  other  two  convection  schemes, 
while  the  second-order  central  difference  scheme  does  not  perform  well  for  either 
steady  or  unsteady  problems  with  strong  convection  effect,  QUICK  is  more 
satisfactory  for  the  unsteady  equation,  yielding  wiggles  mainly  when  the  steady  state 
is  reached.  This  is  because  at  steady  state  the  solution  variation  is  confined  in  a  thin 
boundary  layer  close  to  the  right  hand  boundary,  resulting  in  a  higher  gradient  there. 
Overall,  the  second-order  upwind  scheme  performs  best  if  both  steady  and  unsteady 
solutions  are  needed. 

Figure  2-2  compares  the  solutions  obtained  by  using  the  backward  Euler  and 
Crank-Nicolson  schemes;  they  are  noticeably  different.  For  the  first-order  upwind 
scheme,  with  C  =  1,  the  initial  profile  is  better  preserved  and  the  gradient  is  less 
smeared  out  with    Crank-Nicolson.    In  conjunction  with  the  first-order  upwind 
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scheme,  Crank-Nicolson,  being  less  dissipative  but  yielding  no  wiggles,  performs 
better  than  backward  Euler.  Similarly,  the  performance  of  a  second-order  time 
stepping  scheme  such  as  Crank-Nicolson  is  less  satisfactory  when  combined  with 
other  second-order  convection  schemes  due  to  the  fact  that  both  time  stepping  and 
convection  schemes  are  dispersive.  As  illustrated  in  Figures  2-2  (b)  -(d),  nonphysical 
solutions  are  observed  for  all  three  second-order  convection  schemes.  For  both  the 
central  differencing  and  QUICK  schemes,  the  oscillations  yielded  by  Crank-Nicolson 
are  substantially  larger  than  those  by  backward  Euler.  For  the  second-order  upwind 
scheme,  while  the  Crank-Nicolson  scheme  yields  improved  solutions  with  sharper 
gradients  and  no  oscillations  at  t  =  0.8,  it  produces  a  noticeable  nonphysical 
undershoot  at  t  =  0.3. 

The  results  presented  in  Figures  2-1  and  2-2  are  generic  and  not  affected  by 
the  variations  in  the  number  of  grid  points.  In  order  to  demonstrate  this  fact,  Figure 
2-3  compares  the  performances  between  the  Crank-Nicolson  scheme  and  the 
backward  Euler  scheme,  along  with  the  four  different  convection  approximation 
schemes,  where  the  values  of  P  and  C  remain  the  same  as  before  (10  and  1, 
respectively),  but  the  number  of  nodal  points  is  increased  from  21  to  81.  All  the 
major  characteristics  depicted  in  Figure  2-3  can  be  observed  in  Figure  2-2. 

2.3.2  Nonlinear  Burners  Equation 

For  the  nonlinear  Burgers  equation,  two  different  values  of  viscosity,  v  = 
5x10-'  and  5xl0\  are  used  to  explore  the  effect  of  solution  gradients  on  the 
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performance  of  each  scheme.  Three  time  stepping  schemes,  including  the  first-order 
backward  Euler  scheme,  the  second-order  Crank-Nicolson  scheme,  and  the  third- 
order  Adams-Bashforth/Adams-Moulton  predictor-corrector  scheme,  are  studied 
along  with  the  first-order  and  second-order  upwind  schemes.  Solutions  obtained  at 
two  time  instants,  i.e.,t  =  2  and  3.5, are  presented. 

The  case  of  v  =  5x10"^  is  considered  first.  Figure  2-4  compares  the  solutions 
yielded  by  the  three  time  stepping  schemes  in  conjunction  with  (a)  the  first-order 
upwind  convection  scheme  and  (b)  the  second-order  upwind  convection  scheme.  All 
the  computations  shown  in  Figure  2-4  are  based  on  51  nodes  (Ax  =  0.02),  and  At  = 
0.01,  resulting  in  P  =  4u  and  C  =  0.5u.  Figure  2-4  shows  that  with  first-order 
upwind,  the  numerical  solutions  yielded  by  backward  Euler  and  Crank-Nicolson  are 
fairly  comparable  while  noticeable  improvement  can  be  obtained  by  using  the  third- 
order  A-B/A-M  scheme.  Figure  2-4  indicates  that  a  first-order  scheme  in  space  and 
a  third-order  scheme  in  time  can  actually  form  a  good  combination. 

Regarding  the  use  of  the  second-order  upwind  convection  scheme,  shown  in 
Figure  2-4  (b),  all  three  time  stepping  schemes  perform  well,  producing  solutions  in 
good  agreement  with  the  analytical  one.  In  order  to  measure  more  quantitatively  the 
performance  of  each  individual  scheme  under  study,  Figure  2-5  depicts  an  L-2  error 
norm  with  respect  to  Ax,  where  the  Courant  number  is  held  fixed  for  all  cases,  i.e.. 
At  and  Ax  vary  in  exact  proportion.   The  L-2  error  norm,  CjCt),  is  defined  as  follows, 
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where  </)exact(Ut)  is  the  exact  solution,  ^(i^t)    is  the  numerical  solution,  and  N 

represents  the  total  number  of  nodal  points. 

Several  features  are  revealed  by  Figure  2-5.  First,  with  a  given  convection 
scheme,  the  orders  of  accuracy  of  all  three  numerical  solutions,  with  respect  to  Ax. 
yielded  by  the  three  time  stepping  schemes  are  essentially  the  same.  While  the  local 
truncation  error  analysis  through  a  standard  Taylor  series  expansion  shows  that  the 
numerical  accuracy  of  the  first-order  upwind  convection  scheme  is  proportional  to 
Ax,  Figure  2-5  indicates  that  for  the  present  case  the  errors  are  actually  closer  to 
(Ax)*^  instead  of  Ax.  It  is  noted  that  the  solutions  shown  in  Figure  2-5  are  based  on 
the  fixed  value  of  C,  hence  both  At  and  Ax  contribute  to  the  characteristics  observed 
there;  here  the  overall  accuracy  is  higher  than  that  indicated  by  a  local  truncation 
error  analysis  based  on  the  spatial  accuracy  of  the  first-order  upwind  convection 
alone.  As  to  the  second-order  upwind  convection  scheme.  Figure  2-5  (b)  shows  that 
the  solutions  using  both  Crank-Nicolson  and  A-B/A-M  yield  errors  proportional  to 
(Ax)^  while  that  using  backward  Euler  is  only  slightly  worse.  Hence,  for  the  first- 
order  upwind  convection  scheme,  the  resulting  accuracy  of  all  solutions  are  better 
than  proportional  to  (Ax),  while  for  the  second-order  upwind  convection  scheme,  they 
are  essentially  proportional    to  (Ax)l     Overall,  Figures  2-4  and  2-5  clearly 
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demonstrate  that  the  combined  performance  of  the  discrete  operators  both  in  time 
and  in  space  may  not  be  the  same  as  that  suggested  by  a  local  truncation  error 
analysis  of  the  individual  operator. 

The  case  of  =  5x10"^  represents  a  problem  of  relatively  smooth  solution 
profiles.  A  more  stringent  case  of  v  =  5x10'^  is  considered  next.  Figure  2-6 
compares  the  numerical  solutions  corresponding  to  the  different  combinations  of  the 
convection  and  time  stepping  schemes.  Comparing  Figures  2-4  and  2-6,  it  is  clear 
that  the  degrees  of  satisfaction  yielded  by  these  schemes  between  the  two  values  of 
u  are  very  different.  With  a  reduced  u,  the  solution  gradients  stiffen,  and  the 
numerical  wiggles  are  more  prone  to  appear.  Figure  2-6  (a)  shows  that,  even  with 
the  use  of  the  first-order  upwind  convection  scheme,  which  is  viewed  as  usually  too 
dissipative,  nonphysical  wiggles  develop  when  the  A-B/A-M  time  stepping  scheme 
is  employed.  With  either  backward  Euler  or  Crank-Nicolson,  on  the  other  hand,  the 
solutions  are  of  smooth,  diffusive  characteristics  that  exhibit  no  overshoots.  The 
difficulties  of  resolving  the  sharp  gradients  become  more  pronounced  for  the  second- 
order  upwind  scheme,  shown  in  Figure  2-6  (b),  in  which  the  numerical  wiggles  appear 
in  all  solutions. 

The  L-2  error  norms  of  the  solutions  of  u  =  5x10'*  are  summarized  in  Figure 
2-7  where  the  Courant  number  again  remains  unchanged  as  Ax  is  varied.  Unlike 
those  curves  shown  in  Figure  2-5  for  v  =  5x10'^,  the  results  in  Figure  2-7  for  p  = 
5x10"' are  of  inconsistent  appearances,  depicting  no  recognizable  correlations  with  Ax 
among  different  discretization  schemes  or  between  different  time  instants.  In 
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particular,  the  solutions  associated  with  the  A-B/A-M  scheme  are  most  erratic;  they 
can  either  be  less  accurate  than  the  solutions  using  other  time  stepping  schemes  or 
produce  more  error  as  Ax  is  reduced.  The  results  presented  so  far  clearly  suggest 
that  the  three  time  stepping  schemes  exhibit  different  dispersive  characteristics;  as 
the  leading  local  truncation  error  of  a  time  stepping  scheme  becomes  higher-order 
in  At,  the  resulting  solutions  tend  to  be  more  dispersive,  in  form  of  either  long  wave 
overshoot  or  short  wave  wiggles. 

The  above  discussions  are  based  on  a  single  set  of  the  Courant  and  cell  Peclet 
numbers.  Since  each  of  the  time  stepping  and  convection  schemes  possesses  its  own 
diffusive  and  dispersive  characteristics,  it  appears  useful  to  adjust  the  values  of  P  and 
C  so  that  a  time  stepping  scheme  of  particular  dispersive  and  diffusive  characteristics 
can  be  utilized  in  conjunction   with  a  convection   scheme  of  complementary 
characteristics  to  improve  the  overall  solution  accuracy.  That  such  a  possibility  does 
exist  is  illustrated  in  Figure  2-8,  where  both     5xl0■^and  Ax,  0.02, remain  the  same 
as  those  used  in  Figure  2-6,  while  At  is  increased  from  0.01  to  0.05,  resulting  in  a  C 
=  2.5u.  By  increasing  the  size  of  At,  the  backward  Euler  scheme  exhibits  stronger 
diffusion  characteristics  while  the  other  two  higher-order  time  stepping  schemes 
exhibit  stronger  dispersive  characteristics.  Hence,  as  demonstrated  in  Figure  2-8  the 
first-order  upwind  convection  scheme,  being  generally  more  diffusive,  can  be  used 
along  with  the  Crank-Nicolson  scheme,  a  more  dispersive  scheme.  The  net  outcome 
is  that  they  can  produce  satisfactory  solutions  exhibiting  no  spurious  wiggles  with  a 
larger  time  step,  yielding  saving  in  the  computing  cost.   On  the  other  hand,  the 
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second-order  upwind  convection  scheme,  being  less  diffusive,  can  apparently  benefit 
from  the  backward  Euler  time  stepping  scheme,  which  is  more  diffusive.  These 
interesting  features  suggest  that  the  different  characteristics  in  convection  and  time 
stepping  schemes  can  indeed  be  selectively  combined  to  form  an  overall  algorithm 
capable  of  suppressing  excessive  numerical  dispersion  and/or  diffusion  with 
reasonable  At. 

2.4  Von  Neumann  and  Spectral  Analyses 
In  order  to  gain  some  insight  of  the  dissipative  and  dispersive  characteristics 
of  a  given  numerical  scheme,  the  von  Neumann  stability  analysis  (Warming  and 
Hyett  1974,  Morton  1971)  and  the  discrete  fast  Fourier  transform  (FFT)  spectral 
analysis  are  utilized.  The  von  Neumann  stability  analysis  is  applicable  to  linear, 
constant  coefficient  problems  with  periodic  boundary  conditions. 
The  elementary  solution  for  linear  Burgers  equation  is 

=  g (2.20) 

where  i  =  V-1, 

k„  =  2ir/L„  is  the  wave  number  of  the  m-th  component, 
and  L„  is  the  wavelength  of  the  m-th  component. 

Based  on    the  elementary  solution,  the  amplification  factor  of  the  exact 
solution  can  be  defined 
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G  =  ^(^^^)  =  e'^^'e-'^""  (2-21) 
0(0 


where  |  Gj  |  =  exp(-C/3Vp)  is  the  amplification  magnitude,  0  =  -  jSC  is  the  phase 
angle,  jS  =  k„  Ax,  and  27rAx//3  is  the  wavelength. 

It  is  clear  that  depends  on  the  cell  Peclet  number,  Courant  number,  and 
the  wave  number.  The  amplification  factors  yielded  by  the  von  Neumann  analysis 
for  the  backward  Euler  and  Crank-Nicolson  time  stepping  schemes  along  with  various 
convection  schemes  for  linear  Burgers  equation  are  given  in  Tables  1  and  2. 
respectively. 

In  order  to  illustrate  the  characteristics  of  each  scheme  more  clearly,  a  series 
of  results  are  computed  for  P  =  10^,  C  =  0.5,  Ax  =  0.05,  and  a  periodic  boundary 
condition.  The  initial  condition  is  defined  as 

_/l,  /or  0.05<^  <0.2  or  0.55<;c  <0.7 


'^^'''^^  =  {o,  themlse 


Each  of  Figures  2-9  to  2-16  includes  four  parts  :  (i)  the  comparison  of  exact 
and  numerical  solutions  at  t  =  2At,  (ii)  the  spectral  analysis  of  the  numerical 
solutions  yielded  by  FFT,  and  the  comparison  of  (iii)  the  amplification  factors  as 
well  as  (iv)  the  phase  angles  according  to  the  von  Neumann  analysis.  The 
performance  of  the  backward  Euler  time  stepping  scheme  along  with  different 
convection  schemes  are  presented  in  Figures  2-9  to  2-12;  the  corresponding  results 
of  the  Crank-Nicolson  time  stepping  scheme  are  presented  in  Figures  2-13  to  2-16. 
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For  the  linear  problem  subject  to  a  periodic  boundary  condition,  the  ratios  of 
magnitudes  at  each  wave  number  between  two  consecutive  time  steps  shown  in  the 
spectral  plot  (iii),  are  identical  to  the  amplification  magnitude  shown  in  (iii). 

It  is  noted  that  with  a  high  Peclet  number,  P  =  10\  and  a  periodic  boundary 
condition,  the  exact  solution  essentially  maintains  its  initial  profile  just  like  a  pure 
wave  equation.  Consequently,  the  spectral  distributions  of  the  exact  solution  remains 
unchanged  with  time,and  the  amplification  magnitude  shown  in  (iii)  are  equal  to  1 
at  all  wave  numbers.  Several  observations  can  be  made  by  inspecting  Figures  2-9  to 
2-16.  First,  Crank-Nicolson  does  not  yield  less  dissipation  than  the  backward  Euler 
scheme  for  all  wave  numbers.  In  fact,  except  for  the  case  of  the  central  difference 
convection  scheme,  part  (iii)  of  Figures  2-9  to  2-16  demonstrate  that  for  all  three 
upwind  type  of  convection  schemes,  Crank-Nicolson  is  less  dissipative  than  backward 
Euler  only  in  the  long  wave  length  regime;  with  shorter  wave,  it  is  actually  more 
dissipative  than  backward  Euler.  Secondly,  the  amplification  magnitudes  of  the 
second-order  central  difference  scheme  are  very  close  to  the  exact  values;  when 
combined  with  Crank-Nicolson,  it  is  indistinguishable  from  the  exact  profile,  i.e.  there 
is  no  numerical  damping  created  by  the  scheme,  which  is  good.  However,  with  the 
central  difference  scheme,  the  short  waves  travel  with  much  smaller  phase  speeds 
than  the  exact  solutions.  It  is  these  phase  speed  errors  in  the  short  wavelength  that 
cause  the  central  difference  scheme  to  perform  unsatisfactorily.  This  observation 
explains  why  the  central  difference  scheme  usually  yields  solutions  with  highly 
noticeable  2Ax  oscillations. 
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The  other  interesting  aspect  revealed  by  the  present  analysis  is  that  the 
second-order  upwind  scheme  is  actually  more  dissipative  than  the  first-order  upwind 
scheme  at  short  wavelengths  as  indicated  by  part  (iii)  of  Figures  2-9  and  2-11.  With 
both  Crank-Nicolson  and  backward  Euler,  the  second-order  upwind  scheme  is  less 
dissipative  than  the  first-order  upwind  scheme  only  in  the  long  wave  regime.  This 
feature  may  be  useful  if  a  solution  profile  contains  sharp  gradients  because  the  phase 
speed  errors  associated  with  short  waves  are  less  pronounced.  Furthermore,  by 
comparing  Figures  2-11  and  2-15,  it  becomes  clear  that  it  is  the  phase  errors  in  the 
long  wavelength  that  are  mostly  responsible  for  the  inaccuracies  created  by  the 
combined  Crank-Nicolson/second-order  upwind  schemes.  This  can  explain  why  its 
solutions  exhibit  a  single  long  wave  overshoot  instead  of  2Ax  short  wave  oscillations. 
For  the  first-order  upwind  scheme,  it  appears  that  in  the  long  wave  regime,  Crank- 
Nicolson  is  not  only  less  dissipative  but  also  more  accurate  in  phase  speed. 
Accordingly,  as  already  concluded  previously,  backward  Euler/ second-order  upwind 
or  Crank-Nicolson/first-order  upwind  form  good  combinations.  The  performance  of 
QUICK  lies  between  the  central  difference  and  upwind  schemes. 


2.5  Summary 

Both  linear  and  nonlinear  Burgers  equations  have  been  used  as  the  model 
problems  to  investigate  the  interaction  of  the  time  stepping  and  convection  schemes 
for  unsteady  flow  simulation.  The  von  Neumann  stability  analysis  and  the  FFT 
spectral  analysis  are  also  utilized  to  aid  our  understanding  of  the  performance 
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characteristics  of  these  schemes.  Interesting  features  are  revealed  by  these  analysis. 
For  example,  Crank-Nicolson  is  less  dissipative  than  backward  Euler  only  in  the  long 
wave  regimes;  it  is  actually  more  dissipative  as  the  wavelengths  are  shortened  toward 
2Ax.  Similarly,  with  the  same  time  stepping  scheme,  the  second-order  upwind 
scheme  is  less  dissipative  than  the  first-order  upwind  scheme  mainly  for  the  longer 
waves.  Consequently,  for  the  second-order  upwind  scheme,  backward  Euler  yields 
more  accurate  phase  speed  than  Crank-Nicolson;  while  for  the  first-order  upwind 
scheme,  it  is  Crank-Nicolson  that  is  more  satisfactory  in  phase  speed.  Regarding  the 
second-order  central  difference  scheme,  because  it  contains  very  little  numerical 
damping  for  all  wavelengths,  the  large  phase  speed  errors,  especially  at  short  wave 
lengths,  cause  its  solutions  to  exhibit  pronounced  2Ax  short  wave  oscillations.  In 
contrast,  for  second-order  upwind,  because  of  its  higher  dissipation  rate  at  short 
wavelengths,  it  is  the  phase  errors  of  longer  waves  that  cause  the  noticeable 
inaccuracies;  consequently,  the  solution  error  tends  to  appear  as  a  single  long  wave 
overshoot.  The  performance  of  QUICK  essentially  lies  between  the  central 
differencing  and  upwind  scheme. 

Overall,  among  all  the  schemes  tested,  either  a  combination  of  first-order 
upwind  for  convection  and  Crank-Nicolson  for  time,  or  a  combination  of  second- 
order  upwind  for  convection  and  backward  Euler  for  time  performs  better.  A 
consistent  second-order  accuracy  for  both  time  and  space  discretizations  does  not 
produce  good  results  for  the  time  dependent  flows  with  large  gradients.  The  results 
suggest  that  a  selective  combination  of  the  convection  and  time  stepping  schemes  of 
appropriate  dispersive  and  diffusive  characteristics  can  attain  good  performance. 


35 


a. 

E 
to 


CO 

a. 


a, 

+ 


A. 

w 

I 


a,|(j 


.S 

CO 

0. 


8 

CM 
+ 

+ 

I 


a,|0 


CO. 

.S 

CO 

a. 


cs. 
a,i<N 


CO. 

ST 

CM 

+ 


CM 

a,|cM 

I 

a,7b 
+ 
a. 

CM 

I 


ca. 
.5 


cs. 


CO 

Q,loo 


a,HJ 


a,|tN 


ea. 

a,|oo 

aTTo 
+ 
a. 

CO  I  oo 
+ 

I 


u 
E 

«  J= 

I-  V] 

O 

M  C 

o. 

3 


U 

o 

l_  c 

"2  ^ 
o  t 

■6  ^ 

c 

o  IS 

U  w 
</)  c 

u 


u 

h  E 
u 


u 
■a 


o 

■^^ 

o  .5 
u  > 

00  Q. 


U  E 
—  « 


36 


c 
_o 

u 
E 

CO 


a. 


en. 
u 

ST 

+ 


I 


.5 

V) 

a. 


a, 

+ 

+ 

tN 

+ 

a. 

+ 

CM 


CO. 

a. 


If 


I 
I 

CN 


ca. 
.5 

CO 

a. 


I 

ca. 


a,Tu 
+ 

I 


ca. 
.S 


.S 

09 

a,|cs 

I 


■t- 


a. 

+ 


caL 

CM 


a,|rj 
+ 

aTfu 

CM 

I 

a. 

+ 


.s 

a, 

CN 

I 

ca. 

<N 

a,|<N 


ca. 


a, 

(N 
+ 
«N 

+ 
«N 


a,|cN 
a,7(j 

CN 
+ 

a. 

cn  |cN 

+ 

CN 

>^ 

I 


C5 


a,|cN 
+ 

CN 


ca. 


a.|co 

■  + 

a-Tu 

CN 

I 

a, 

cn  |oo 
+ 

CN 


ca. 
.S 

CO 


ca. 

CN 

.s' 

a,|co 


u 

aTTcN 
+ 


ea. 

CN 

8 
a,|oo 


a,|0 

CN 
+ 

a, 

cn  I  oo 
+ 

CN 

I 


II 


E 
u 


E 

u.  U 
U  f 

■a  u 

k.  (/] 
o 

CO  C 

£  S 
Q. 

3 


U 

u.  C 


T3 
C 

u 
(J 


u 

II 

r  « 


o  .5 

M  D. 

3 


U  E 
—  « 


37' 


5 

s 

II 


d 


in 


o 


■S  °  "  T  « 
O  01  o  <-> 

-  S  •-  u  ^ 
=  i3)  -o  tn  a 

CM 

d 


o 


^  - 

«  II 

u  5 


c 
o 


CO 


3  Wl 


00 

d 

(A 

O 

u 
j= 
o 

II 

c 

ISO 

on 

C 

'5 

w 

f8 

u 

Q. 

c 

u 

E 

> 

"(3 

c 

•o 

o 
U 

CO 

no 

CM 

a 


s 
on 


38 


—  Backward  Euler 

-  -  Crank-Nicolson 


V 

E 

u 


c 

'o 
c 
u 


■a 

Urn 

c 

u 
u 


u 
•13 


13 

c 
o 
u 
u 


u 

E 
u 

u 


0. 

3 


13 


1« 


a 
II 

( 

r>  / 
0  / 

1 

T3 

00 
d 


u 

E 
u 

u 

y 
5 


u 
c 

Q. 

3 


1) 
"O 

o 

I 

•o 
c 
o 
u 
u 

CO 


00 

e 


•8 


O  4J 
§  J 

—  o 

5*  00 

«  c 

en  & 
^  Q. 

C 

ca  *" 

u 

=  1 

>»  o 
■o  o 

<u  Z 


C/5 


C 


i  = 
.2 

1 1 

I- 

O  V 


3 

as 


39 


41 


43 


44 


4) 


U 


> 
C 

o 
u 


5 

a. 


c 
o 
u 
u 


u 
E 


u 

> 
c 
o 
o 


a. 


o 


o 


O  O 

^  II 

L- 

5  < 

u 

CO  "O 

c 


o 
d 

II 

X 

<3 


c 
o 


13  w-i 


OA  ^ 

•|:£ 


u  o 
5 

uu  aa 


oo 

t 

u 

w 

3 
00 


45 


46 


47 


d  '■^      o      d      o  a 

IP 


50 


52 


<u 

V. 
CO 


LO 

d 


2  if;  o  in  o  ui  o"* 

-.         .  (M     ^     X  Q 

o  o  d  d  d  d  d  d 


■i.  =; 
2  o  *  < 

il  M 
"    "    K  I. 

U  <  _ 


in 
d 


re 

E 


CO 

u 


E 

3 
00 


s  "I 

c 

CO 

3 
J  II 

8  a, 
z 

5 

=  "i: 
. — 

^  § 
u 

E  >. 
■  — 

^  CO 

il 

«  o 

u  ^ 

<r  « 

o  -o 
«  .2 
.Si  S 

oj  — 

CO 


C/5 


£  II 

a.  Q  < 


CHAPTER  3 

MENISCUS  CHARACTERISTICS    WITH  APPLICATION  TO 

THE  EFG  PROCESS 


3.1  Introduction 

As  already  described,  during  the  course  of  fibre  growth,  a  meniscus  is  formed 
adjacent  to  the  solidified  material.  In  order  to  obtain  crystals  of  definite  diameters 
with  given  growth  rates  or  rates  of  pulling,  the  control  of  the  behavior  of  the 
meniscus  of  liquid  bridging  the  elevated  solid  and  the  melt  plays  a  key  role.  In  this 
chapter,  we  intend  to  examine  the  behavior  of  the  meniscus  and  shed  some  light  on 
fundamental  aspects  of  generic  interest  that  have  hitherto  received  only  diffuse 
attention. 

Studies  of  meniscus  profiles  and  properties  have  been  made  by  various 
investigators.  In  the  works  of  Katoh  et  al.  (1990),  and  Pitts  (1976),  and  Huh  and 
Scriven  (1960),  the  meniscus  is  formed  by  dipping  solids  of  different  geometries  into 
baths  of  infmite  extent.  On  the  other  hand,  in  the  crystal-pulling  literature, 
theoretical  treatments  of  meniscus  shape  and  stability  abound  (Michel  1981,Surek 
et  al.  1980,  Tatarchenko  and  Brenner  1980),  mainly  in  the  form  of  linear  analysis, 
and  simulations  of  crystal  growth  processes  including  heat  transfer  have  been  made 
(Backman  1992,  Thomas  et  al.  1986),  with  the  consideration  of  equilibrium  meniscus 
shape.      However,   the   physics  involved   in  this   process,   in  particular,  the 
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thermodynamics  governing  the  menisci  and  the  issues  related  to  contact  angles  have 
not  been  adequately  clarified  (Dussan  V.  1979).  Furthermore,  to  our  knowledge,  the 
possibility  of  multiple  meniscus  shapes  for  a  single  contact  angle  and  consequences 
thereof  has  not  been  studied.  It  is  generally  accepted  (although  not  without 
argument)  that  both  from  considerations  of  force  and  energy,  the  equilibrium  contact 
angle  may  be  obtained  from  the  well  known  Young's  formula. 

Under  debate  is  the  concept  of  the  dynamic  contact  angle,  and  that  of 
advancing  and  receding  contact  angles  (Miller  1985,  Ngan  and  Dussan  V.  1982, 
Schwartz  and  Tejada  1982).  Dynamic  contact  angles  are  a  feature  of  concern  since 
in  such  processes  the  solid  to  which  the  meniscus  is  connected  is  subjected  to  pulling 
as  well  as  rotation  while  being  retracted  from  the  melt.  Specifically,  one  seeks  to 
determine  a  functional  form  for  the  variation  of  contact  angle  0,  with  the  velocity  of 
the  contact  line.  While  success  has  eluded  theoretical  treatment,  experiments  have 
provided  a  broad  picture  of  the  behavior  of  the  apparent  contact  angles  with  contact 
line  motion.  Hysteresis  of  dynamic  contact  angles  implies  that  the  behavior  of  (f>^ 
depends  on  the  direction  of  motion  of  the  contact  line.  The  static  contact  angle 
achieved  at  a  contact  line  that  resulted  from  advancing  fluid,  say  4>^  is  greater  than 
that  formed  by  receding  liquid,  say  (f)^.  Thus,  the  apparent  angle  assumed  under 
static  condition  by  the  liquid-gas  surface  at  the  trij unction  is  either  or  <^r  . 
depending  on  how  it  was  established.  When  the  contact  line  is  compelled  to  move, 
it  will  remain  pinned  while  the  contact  angle  (  via  a  series  of  jumps)  traverses  the 
gap       to  <^R  or  vice-versa,  depending  on  the  direction  of  impending  motion. 
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Thereupon,  the  contact  line  will  move,  and  a  slip  condition  (de  Gennes  1991,Koplik 
et  al.  1989)  and  caterpillar  motion  are  believed  to  be  the  characteristics  that 
accompany  the  translation  of  the  contact  line.  It  is  generally  believed  that  the 
apparent  contact  angle  monotonically  increases  above  <Aa  for  advancing  contact  lines 
and  decreases  below  4>f^  for  receding  contact  lines,  although  the  physics  of  the 
motion  or  an  appropriate  model  (Joanny  and  de  Gennes  1984,  Jansons  1986)  is  not 
established. 

In  the  following  we  will  draw  attention  to  some  of  the  above  issues  concerning 
the  meniscus  formation  caused  by  a  pulling  process.  Although  satisfactory  solutions 
to  these  theoretical  weak  spots  are  awaited,  an  attempt  is  made  to  cast  them  in 
perspective.  Herein,  we  once  again  solve  the  Laplace- Young  equation  for  the 
equilibrium  meniscus  profiles.  The  boundary  conditions  are  chosen  to  correspond 
to  that  of  the  edge-defined  film-fed  growth  (EFG)  process.  We  seek  to  identify  the 
forms  and  properties  of  menisci  and  the  parameters  affecting  shape  and  stability.  In 
particular,  the  Bond  number,  contact  condition,  aspect  ratio  (meniscus  height/base 
radius)  and  the  degree  of  pressurization  (interacting  with  gravity)  significantly 
influence  meniscus  behavior.  Our  motive  is  to  understand  the  equilibrium  behavior 
of  menisci  from  a  fundamental  viewpoint.  As  will  be  demonstrated,  the  issue  of 
equilibrium  meniscus  behavior,  even  without  consideration  of  dynamic  contact  angles, 
holds  many  points  in  need  of  elucidation. 
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3.2  Configuration.  Assumptions  and  Formulation 
The  configuration  chosen  here  is  shown  in  Figure  3-1,  where  h„  re,  and  rt  are 
the  meniscus  height,  radius  of  wetted  region  at  top  surface,  and  die  radius 
respectively.  We  are  particularly  interested  in  the  axisymmetric  situation  related  to 
forming  thin  fibers  (a  few  mm  or  smaller  in  diameter).  As  will  be  discussed,  the 
diameter  of  the  fibre  can  be  quantified  in  relation  to  a  characteristic  capillary  length 
scale  of  a  given  material  through  the  Bond  number.  Since  the  geometry  is 
axisymmetric,  the  appropriate  Laplace- Young  equation  for  the  liquid-gas  interface 
is  applied,  which  in  dimensional  form  is  (Surek  et  al.  1980,  Pitts  1976,  and  Huh  and 
Scriven  1960) 


=  Ap5z  -  Ap 

where  f  is  the  radial  co-ordinate  defining  the  surface  of  the  meniscus,       is  the 
surface  tension  of  the  liquid-gas  interface,  Ap  is  the  difference  in  density  between 
liquid  and  gas  phases,  and  Ap  is  the  difference  in  applied  pressure  between  the 
liquid,  pi  ,  and  gas  phases,  Pg  ,  i.e.,Ap  =  Pi  -  Pg  .  In  our  results,  below,  we  will 
express  Ap  as  a  percentage  of  ambient  pressure,  considered  fixed  at  1  atmosphere. 
Consider  a  capillary  length  defined  as 
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a=  JJl.  (3.2) 

where,  for  example,  for  sapphire  (AI2O3),  7,g  =  0.69  N/m,  Ap  =  3.4x10?  kg/m^  ,  and 
hence  the  capillary  length  is  about  4x10"^  m.  The  Bond  number,  which  relates  the 
radius  of  the  base  or  die  to  the  capillary  length,  and  controls  the  relative  importance 
of  gravity  to  capillary  forces  is  given  by  ,        ^  - 


Bo=^  (3-3) 
a' 


We  specify  the  aspect  ratio,  hjTy,  and  Bond  number,  Bo,  in  each  case.  The  meniscus 
height  is  calculated  from  the  given  aspect  ratio  and  r^.  We  non-dimensionalize  all 
physical  quantities  as  follows,  =^'s representing  non-dimensional  values,  z*=  z/a  ,r  = 
f/a  ,  Ap*=  Ap/(Apga),  Ap  '  =  1.  Again,  using  the  sapphire  as  an  example,  the 
reference  pressure  p,  =  Ap  g  a  is  about  130  N/m^  . 

Thus,  the  non-dimensionalized  profile  equation,  after  having  dropped  all  *'s 

is, 

fzz  1  . 

 T  r  =  ^-^p  (3.4) 


This  equation  is  seemingly  invariant  with  respect  to  Bond  number,  since  the 
Bond  number  does  not  explicitly  appear  in  Eq.  (3.4).  This  is  a  misleading  impression, 
since  the  Bond  number  does  influence  meniscus  profiles  via  the  term  z,  which  is  the 
non-dimensionalized  gravity  term.  For  a  given  aspect  ratio,  h,/rb,  the  effect  of  gravity 
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increases  with  Bo,  since  a  higher  Bond  number  implies  a  taller  meniscus.  When 
Ap?iO,  the  effect  of  Bond  number  will  influence  menisci  through  the  relative  values 
of  hydrostatic  pressure,  Apgz,  and  externally  imposed  pressure  difference,  Ap,  as  will 
be  shown  from  our  calculations. 

The  Young's  contact  condition  at  the  top  of  the  meniscus  is, 

where  7,g,  7,,  and  7,5  are  respectively  the  surface  tensions  of  the  solid-gas,  solid-liquid 
and  liquid-gas  interface.  is  the  static  contact  angle  measured  with  respect  to  the 
liquid-solid  interface. 

It  may  be  noted  that  this  relation  is  strictly  applicable  to  a  perfectly  smooth 
surface  without  energy  heterogeneity.  It  merely  defines  the  static  contact  condition 
at  the  three-phase  line  for  an  ideal  surface.  However,  if  the  solid  surface  is  taken 
to  consist  of  'patches'  of  ideal  smooth  surfaces,  we  may  impose  local  Young's  contact 
angle  at  the  contact  line.  This  situation  has  been  shown  previously  to  offer  one 
plausible  model  for  hysteresis  (Huh  and  Mason  1977,  Eick  et  al.  1975,  Neumann  and 
Good  1972).  It  is  also  evident  that  such  heterogeneities  or  those  associated  with 
crystalline  orientations  or  defects,  may  cause  the  meniscus  to  deviate  from 
axisymmetry.  Contortions  of  the  contact  line  are  developed  in  order  to  ensure  energy 
minimization,  by  locally  enforcing  on  the  liquid-gas  surface  the  Young's  angle  at  the 
contact  line.  This  in  turn  distorts  the  meniscus  surface  which,  for  equilibrium,  has 
to  conform  locally  to  the  Young's  condition.  In  our  work,  we  avoid  the  complications 
associated  with  these  aspects  by  assuming    an  ideal  surface  at  the  top.    It  is  a 
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conventional  practice  to  account  for  the  contact  angle  hysteresis  by  adopting  a 
relation  such  as  Wenzel's,  namely 

where  rf  is  a  roughness  factor  and  (i>,„  is  the  apparent  contact  angle  (Oliver  et  al. 
1977).  So  as  not  to  dilute  our  focus  on  the  fundamental  meniscus  behavior  under 
equilibrium  conditions,  this  practice  is  not  adopted.  It  is  possible,  however,  by 
scanning  through  the  available  range  of  4>o  for  solutions  to  the  Laplace- Young 
equation,  to  incorporate  the  effect  of  Eq.  (3.6)  by  resorting  to  the  appropriate  value 
of  4>,^.  Hence,  a  dynamic  process  with  varying  (^,pp  can  be  treated  simply  as  a 
juxtaposition  of  equilibrium  solutions  with  varying  0<,- 

The  Laplace- Young  equation  being  an  ordinary  differential  equation  of  second 
order,  two  boundary  conditions  are  necessary.  For  the  EFG  process,  the  radius  is 
fixed  at  the  base.  The  liquid-vapor  line  is  assumed  to  be  pinned  at  the  edge  of  the 
die  as  shown  (Surek  et  al.  1980,  Oliver  et  al.  1977),  such  pinning  having  been 

observed  in  practice. 

In  the  present  simplified  framework,  i.e., without  explicitly  accounting  for  heat 
flow  and  solidification,  we  choose  to  picture  the  meniscus  as  attached  to  a  flat  plate. 
Thus,  at  the  top,  for  a  given  meniscus  height,  h,  ,  the  contact  angle  is  specified  and 
the  radius  of  the  trij unction  point  can  be  computed.  This  is  in  cognizance  of  the  fact 
that  the  Young's  contact  angle  is  a  material  property,  independent  of  the  geometry 
of  the  meniscus.  Such  a  boundary  treatment  also  offers  a  simple  way  to  simulate 
fibre  diameter  variations.  However,  the  fact  that  a  Neumann  boundary  condition  is 
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specified  at  the  top  results  in  the  existence  of  multiple  solutions.  In  fact,  it  is  known 
that  a  boundary  value  problem  generally  does  not  guarantee  a  mathematically  unique 
solution.  Hence,  guidance  will  be  needed  to  help  select  the  solution  that  is  physically 
stable.  The  implications  of  mathematically  non-unique  solutions  on  the  realizability 
of  the  meniscus  shape  will  be  discussed  later. 

3.3  Energy  Formulation 
The  thermodynamics  of  the  meniscus  may  be  studied  by  defining  the  energy 
change  in  forming  the  meniscus.  The  total  energy  E  of  this  system  is  given  by  the 
sum  of  the  gravitational  potential  energy  and  the  free  energy  of  the  surface.  To 
nondimensionalize  energy, 

£•  =  — ^  (3.7) 
^pga* 

Total  energy  resulting  from  the  formation  of  the  meniscus,  E,  is  equal  to  the 
potential  energy  (from  the  effective  head)  E,  plus  surface  energy  of  meniscus  E,  plus 
surface  energy  of  the  wetted  area  of  the  top  plate  E,.  The  first  component  of  the 
energy,  that  due  to  gravity  and  imposed  pressure  is, 

£i  =  fitf^  (Apgz+A/?)  ek  (3-8) 
The  second  component,  corresponding  to  that  expended  in  forming  the  free  surface 
is, 
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(3.9) 


The  third  component,  the  energy  of  wetting  of  top  surface  is, 


^3  =  (Yfa- Y,^)  T^r] 


(3.10) 


where  r,  is  the  radial  location  of  the  trijunction  at  the  top  of  the  meniscus. 

We  follow  the  usual  practice  of  replacing  the  difference,  (7,,  -7,5)  in  the  third 
component  using  the  Young's  contact  condition  (Pitts  1976,  and  Huh  and  Scriven 
1960).  This  procedure  facilitates  evaluation  of  this  component,  since  the  values  of 
contact  angle   4>q  and       are  often  the  only  experimentally  obtainable  quantities. 
Thus,  the  third  component  of  energy  may  be  written  as 


The  total  energy  then  is,  after  nondimensionalizing  and  dropping  all  *'s. 


By  scanning  through  the  range  of  possible  solutions  the  total  energy  E  can  be 
computed  as  a  function  of  0„  the  contact  angle  which  is  contained  in  the  solution  for 
meniscus  shape,  including  that  corresponding  to  the  specified  value  4)^.  The  decision 
regarding  stability  of  the  meniscus  can  then  be  arrived  at  based  on  whether  the 
specified  value  of  4>o  minimizes  E  or  not. 


£3  =  Y/^  cos4)^  sir] 


(3.11) 


(3.12) 


+  j27iy(lH//)'/2^  +  7ir;cos<t), 
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3.4  Numerics 

The  Differential  equation  governing  the  meniscus  shape,  the  Laplace- Young 
equation,  is  nonlinear  and  has  to  be  solved  numerically.  The  integration  of  the 
nonlinear  ordinary  differential  Eq.  (3.4)  is  performed  using  the  fourth-order  Adams 
predictor-corrector  method,  using  the  Runge-Kutta  method  for  starting.  One- 
hundred  and  one  points  were  used  along  the  z-  direction,  as  deemed  sufficient  from 
grid  retlnement  studies.  The  procedure  adopted  was  to  use  the  shooting  method 
beginning  from  the  fixed  base  at  the  edge  of  the  die.  Shooting  was  performed  at 
various  starting  angles,  henceforth  0b  •  The  resulting  slopes  at  the  top  of  the 
meniscus  were  then  stored.  For  any  specified  slope  at  the  top,  as  required  for  a 
given  contact  angle,  the  profile  was  determined  by  using  a  bisection  method  by 
sweeping  through  the  solutions  stored  from  the  initial  value  treatment.  The  integrals 
in  the  expression  for  energy  were  computed  using  the  trapezoidal  rule. 

3.5  Results  and  Discussion 
As  previously  mentioned,  the  meniscus  profiles  are  significantly  influenced  by 
Bond  number  (Bo),  applied  pressure  (Ap),  aspect  ratio  (h^  /T^,),  and  the  value  of 
contact  angle.  Our  procedure  was  to  specify  the  Bond  number  and  aspect  ratio, 
thereby  obtaining  the  value  of  rb  from  Eq.  (3.3)  ,  and  meniscus  height  from  aspect 
ratio  and  rb. 

As  already  mentioned,  an  interesting  feature  of  our  solutions  for  this  model 
of  Edge-defined  Film-fed  growth  is  the  presence  of  multiple  solutions  for  a  given 
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contact  angle.  Unless  a  methodical  sorting  is  performed,  a  numerical  solution 
procedure  may  not  pick  up  all  the  possible  solutions.  In  the  work  of  Katoh  et  al. 
(1990),  multiple  solutions  seem  to  have  been  obtained  only  in  certain  cases.  This 
may  be  due  to  their  operating  regime  of  parameters,  such  as  Bond  number,  aspect 
ratio  etc.  From  our  solutions  it  appears  that  the  tendency  to  obtain  multiple 
solutions  and  interesting  profile  shapes  is  greater  for  taller  menisci,  and  is  also 
influenced  by  the  applied  pressure. 

3.5.1  Basic  mMeniscus  Behavior  with  Upward  Pulling 

Consider  a  typical  case  (Case  1)  shown  in  Figure  3-2.  The  parameters  are  Bo 
=  4x10  "^  h,/rb  (aspect  ratio)  =  0.5,  Ap  =  0.  The  curves  corresponding  to  0^ 
against  4)^,,  and  r,/rb  against  (t>^  clearly  show  the  multiple  solution  behavior.  Thus, 
there  is  only  an  interval  of  possible  contact  angles  0c  that  can  be  achieved  at  the 
top  of  the  meniscus.  Outside  this  range  of  angles  no  solutions  exist.  This  is  not  due 
to  instability  of  possible  equilibrium  profiles,  but  because  it  is  not  even  possible  to 
obtain  profiles  satisfying  the  Laplace- Young  equation.  Also  shown  in  Figure  3-2  are 
sample  meniscus  profiles,  corresponding  to  (f>o  of  60°  and  90°,  respectively. 

As  shown  in  Figure  3-2,  when  the  value  of  0o  is  specified,  and  E  is  calculated, 
distinct  extrema  are  obtained  at  specified  values  of  <t>o.  Among  the  multiple 
solutions,  the  one  corresponding  to  a  minimum  on  the  curve  is  the  physically 
realizable  stable  meniscus  profile.  When  either  a  maximum  or  a  non-extremum  point 
arises,  we  may  surmise  that  such  solutions  belong  to  the  unstable  branch. 
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That   the   Young's  condition    is  in  fact  the  condition    for  stability  is 
straightforward  to  show,  say  in  the  case  of  a  sessile  drop.  In  such  a  (constant  volume) 
case,  when  gravity  is  ignored,  the  drop  assumes  the  shape  of  the  arc  of  a  circle  (Dyson 
1978).    Then  the  minimization   result  becomes   merely  a  geometric  condition 
compatible  with  minimum  surface  area.   When  gravity  is  considered  the  drop  no 
longer  remains  the  element  of  an  arc  ,  since  minimization  of  surface  free  energy  is 
not  sufficient,  it  being  necessary  to  minimize  the  total  energy  including  gravitational 
potential.    In  the  case  of  menisci  above  infinite  baths,  Pitts  (1976)  and  others  have 
shown  the  minimization  principle  via  variational  methods.    It  is  the  conventionally 
held  view  that  Young's  contact  condition  can  be  obtained  both  from  consideration  of 
forces  as  well  as  energy  (however,  see  Johnson  (1959)  and  discussion  therein  for 
criticism  of  this  view).    However,  it  is  interesting  to  see  that  in  the  present 
configuration,  multiple  shapes  corresponding  to  energy  minimum,  maximum  and  even 
non-extremum  may  all  be  observed  at  the  given  contact  angle  4>o.  Furthermore, 
solutions  to  the  Laplace- Young    equation,  which  is  after  all  an  equilibrium  relation 
across  the  interface,  are  all  equilibrium  solutions.    From  the  standpoint  of  local 
energy  minimization,  there  may  be  more  than  one  stable  shape.    Thus  the  end 
condition  and  the  path  of  development  of  the  meniscus  selects  the  stable  solution. 
However,  in  a  process  in  which  the  aspect  ratio  is  varied,  an  unstable  meniscus  may 
form  first.  This  can  happen  if  the  new  trijunction  location  corresponding  to  this  shape 
is  closer  to  the  initial  trijunction  location.    Then,  in  achieving  the  stable  solution 
profile  for  the  given  height,  one  may  pass  through  an  unstable  solution  profile. 
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As  discussed  in  connection  with  the  non-dimensional  form  of  the  Laplace- 
Young  equation,  for  low  Bond  numbers,  the  term  corresponding  to  gravity  is  small 
compared  to  the  curvature  terms  on  the  left  hand  side  of  Eq.  (3.1).  Thus  the  Bond 
number  serves  principally  to  relate  the  dimensions  of  the  meniscus  to  the  capillary 
length  scale  according  to  Eq.  (3.3).  However,  for  higher  Bond  numbers,  the 
hydrostatic  term  on  the  right  hand  side  of  Eq.  (3.1)  becomes  comparable  to  the 
curvature  terms  on  the  left  hand  side,  and  gravity  does  significantly  impress  on 
meniscus  profiles,  as  will  be  illustrated  in  the  results  to  follow. 

For  taller  menisci,  i.e.,  higher  aspect  ratios,  solutions  to  menisci  appear  to 
offer  a  wider  variety  of  shapes  and  range  of  angles  as  well  as  a  greater  tendency  to 
multiple  solutions.  Retaining  Bo =4x10'^  and  Ap=0,  as  in  Case  1 ,  we  can  increase 
the  height  of  the  meniscus.  For  example,  for  hjr^,  =2,  which  represents  a  taller 
meniscus,  the  available  range  of  contact  angles  will  shrink.  Thus,  as  expected,  when 
the  height  of  the  meniscus  is  increased,  in  the  absence  of  pressurization,  a  large 
number  of  solutions  to  the  Laplace- Young  equation  fail  to  exist. 

Cases  2  and  3  (Figs.  3-3  and  3-4  respectively)  are  illustrative  of  the  effects  of 
pressure  and  Bond  number  on  menisci  with  a  higher  aspect  ratio  where  solutions 
with  selected  contact  angles  are  shown.  For  Case  2,  Bo=2.5xlO'*,Ap=5%,  hjr^  =2. 
For  this  low  Bond  number  case,  from  Figure  3-3,  we  notice  that  for  a  contact  angle 
of  16°,  three  equilibrium  profiles  are  obtained.  The  possible  contact  angles  at  the 
top  were  in  the  range  of  18°  to  26°.  These  profiles  may  not  be  seen  in  a  physical 
situation,  since  the  contact  angles  corresponding  to  the  materials  of  the  system  may 


66 

not  lie  in  this  range  at  all.  Thus,  ultimately  the  existence  of  menisci  will  be  dictated 
by  the  range  of  contact  angles  <l>^  at  the  top  permitted  by  the  materials,  the  range  of 
0b  for  which  the  pinning  condition  is  maintained  at  the  edge  of  the  die,  and  the 
stability  of  the  equilibrium  solution  (of  the  Laplace- Young  equation)  obtained. 

The  contrasting  situation  to  Case  2  emerges  in  Case  3  (Figure  3-4).  Here  the 
Bo  =  2.5x10^ (h,/rb  =  2,and  Ap  =  5%  as  in  Case  2)  is  10^  times  that  of  Case  2.  The 
available  range  of  contact  angles  is  considerably  widened,  and  a  greater  variety  of 
meniscus  profiles  than  for  the  lower  Bo  in  Case  2  can  be  obtained.  Comparison  of 
profiles  with  those  in  Figure  3-3  shows  the  difference  that  Bond  number  makes  to 
the  effect  of  pressure.  Changes  in  curvature  arise  in  the  profile.  This  implies  that 
for  the  given  parameters,  on  the  right  hand  side  of  Eq.  (3.1),  terms  of  pressure  and 
gravity  must  compete  to  determine  the  shape  of  the  meniscus.  Thus,  as  mentioned 
before,  for  Ap  0,  the  Bond  number  can  significantly  influence  meniscus  shapes  and 
permissible  solutions.  This  can  be  seen  from  the  scaling  that  we  have  employed  in 
the  non-dimensionalization.  With  fixed  h,  /rt,  Eq.  (3.3)  shows  that  the  height  of  the 
meniscus  is  scaled  with  respect  to  Bond  number  as  h,  «  y/Bo.  Hence,  for  given  Ap, 
g  and  a,  the  meniscus  height  increases  with  Bond  number.  For  higher  Bo,  by 
definition,  the  relative  importance  of  gravity  (Apgy)  becomes  more  significant 
compared  with  the  pressure  term.  Thus,  the  competition  between  gravity  and 
pressure  in  the  force  balance  leads  to  a  variety  of  curvatures  for  higher  Bond 
numbers.  For  lower  Bo,  a  relatively  small  Ap  can  overwhelm  the  gravity  term,  and 
the  curvature  is  less  sensitive  to  location  along  the  z-axis. 
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Another  aspect  to  note  in  Figs.  3-3  and  3-4  is  that  in  the  E  vs.  <t>^  curves  both 
maxima  and  minima  are  obtained  at  specified  <t)„.  As  a  contrast,  however,  consider 
Case  4  (Figure  3-5).  Here,  we  have  Bo=  4x10^,  Ap=5%  and  with  the  h,  /r^, 
increased  to  5.  While  the  other  properties  appear  to  follow  the  general  trend  of 
Case  3,  the  behavior  of  energy  is  distinctive.  Here,  for  (t)„  =  90°,  both  minimum 
and  maximum  are  obtained.  However,  the  minimum  energy  state  so  obtained  is  only 
local.  There  appears  to  be  a  global  minimum  at  130°.  For  4>„  =  130°,  Figures  3-5 
(c)  and  (d),  the  local  minimum  occurs  at  around  100°  and  there  appears  to  be  a  non- 
extremum  point  corresponding  to  the  specified  indicating  that  although  only  one 
meniscus  profile  exists,  it  is  not  a  stable  one.  Thus,  for  this  case  the  decision 
regarding  the  stable  meniscus  profile  becomes  more  involved. 

In  summary,  the  numerical  simulation  reveals  that  depending  on  the  range  of 
parameters  a  variety  of  solutions  are  possible  for  a  given  contact  angle.  Multiple 
energy  minima,  implying  multiple  stable  shapes,  or  no  minima,  implying  unstable 
shapes  may  be  obtained  for  a  given  contact  angle. 

3.5.2 Effect  of  Pull  Direction  on  Meniscus  Behavior 

The  issues  of  the  selection  of  a  stable  meniscus  shape  out  of  the  possible 
solutions,  and  the  implications  of  variation  of  h,/rb,  as  well  as  the  direction  of  pulling 
on  meniscus  stability,  can  be  better  addressed  through  the  energy  considerations 
stated  earlier.  In  studying  the  variation  of  meniscus  shape  and  top  radius  with  height 
the  effect  of  the  following  factors  is  investigated.  They  are  upward/downward  pulling 
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with  respect  to  direction  of  gravity,  Bond  number  (Bo),  pressure  difference  between 
the  liquid  and  gas  phases  (Ap),  and  aspect  ratio  (h,/rb). 

The  results  to  be  discussed  in  the  following  are  obtained  by  fixing  rb  and 
varying  hjr^  from  0.5  to  2.5  upward  and  downward,  respectively.  Solutions  to  the 
Laplace- Young  equation  are  obtained  for  each  value  of  hjTf,  and  the  range  of 
available  solutions  in  terms  of  (^<.  and  resulting  radius  at  the  top  (expressed  as  r,/rb) 
are  determined.  Figure  3-6  compares  the  characteristics  of  the  meniscus  for  the 
upward  and  downward  pulling  situations,  with  Bo  =  10'^  and  Ap  =  0.  The  figures 
display  the  solution  range  of  <t>^  versus  h<./rb,  the  solution  range  of  r^rb  versus  h,/rb, 
and  possible  meniscus  profiles  for  a  prescribed  contact  angle  <Ao  =  10°.  As  shown, 
in  the  range  of  0^  for  which  the  Laplace- Young  equation  does  offer  solutions,  the 
possible  maximum  contact  angle  (f),  decreases  with  increasing  hjr^,,  while  the  possible 
minimum  contact  angle  (t)^  changes  little  during  the  pulling.  Consequently,  the 
solution  range  of  shrinks  when  h,/rb  is  increased,  implying  that  unless  the 
materials  comprising  the  system  are  such  that  the  contact  angle  at  the  top  has  a  low 
value,  the  meniscus  detaches  from  the  solid  surface  when  hjvy,  is  increased. 

It  is  interesting  to  note  that  with  no  pressurization,  while  the  direction  of 
pulling  causes  little  difference  in  permissible  range  of  contact  angles,  it  affects  the 
trij  unction  location  in  a  more  profound  manner.  The  range  of  permissible  trij  unction 
locations  is  much  narrower  for  upward  pulling  than  for  downward  pulling.  This 
feature  is  illustrated  by  the  selected  meniscus  profiles  corresponding  to  4>^  =  10°, 
shown  in  the  bottom  of  Figure  3-6. 
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With  Bo  =  10"^  when  hjr^  is  less  than  1,  the  hydrostatic  pressure  is  small 
relative  to  the  curvature  terms  on  the  right  hand  side  of  the  Laplace- Young  equation. 
Hence  the  meniscus  profiles  are  essentially  determined  from  the  geometrical 
constraints  alone,  meaning  that  within  this  regime  the  specification  of  r^,  at  the 
bottom  and  <i>^  at  the  top  dictates  the  meniscus  shape.  Figure  3-6  clearly  indicates 
that  the  effect  of  pulling  direction  is  noticeable  only  when  the  aspect  ratio  is  high 
enough.  The  critical  value  of  h./rt  for  the  hydrostatic  pressure  to  exert  a  significant 
effect  is  Bond  number  dependent.  It  appears  that  the  more  appropriate  indicator  of 
the  importance  of  hydrostatic  pressure  is  (h^/rJ^Bo,  i.e., an  effective  Bond  number 
defined  according  to  h,  instead  of  rb. 

We  next  focus  on  the  interaction  of  the  direction  of  pulling,  level  of 
pressurization  and  variation  of  h,/rb.  Figure  3-7  depicts  the  solution  characteristics 
for  the  case  of  Bo  =  10■^  Ap  =  1%,  i.e.,  there  is  an  excess  pressure  inside  the 
meniscus.  Regarding  the  overall  range  of  permissible  meniscus  solutions,  both  in 
terms  of  0^  and  r,,  the  direction  of  pulling  does  not  make  any  noticeable  difference, 
seemingly  indicating  that  the  externally  imposed  pressurization  dominates  the 
hydrostatic  pressure.  However,  depending  on  the  value  of  h,/rb  the  detailed  meniscus 
shapes  can  be  significantly  affected  by  the  pulling  direction.  For  example,  with  h,/rb 
less  than  or  equal  to  1,  little  difference  is  found  between  meniscus  shapes  formed  by 
upward  and  downward  pulling.  For  h,/rb  of  2  or  higher,  however,  the  differences  can 
be  substantial.  Specifically,  for  h,/rb  =  2,  with  upward  pulling,  three  possible 
meniscus  shapes  exist,  while  with  downward  pulling,  only  one  profile  can  be  obtained. 
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When  h./fb  is  raised  to  2.5, for  either  direction  of  pulling  only  one  meniscus  profile 
is  obtained,  but  the  location  of  trijunction  point  for  the  downward  pulling  case  differs 
from  that  of  the  upward  pulling  case. 

The  energy  profiles  corresponding  to  the  cases  in  Figure  3-7  are  presented  in 
Figure  3-8;  for  all  the  values  of  h,  /rt  considered  here,  both  directions  of  pulling  yield 
energy  profiles  that  appear  essentially  identical.  The  energy  curve  for  h^  /t^,  =  2 
shows  that  for  the  upward  pulling  cases,  all  three  possible  meniscus  shapes  contain 
virtually  the  same  amount  of  energy  despite  their  noticeable  differences  in  trijunction 
location.  Significantly,  with  hJVi,  of  2  or  higher,  the  energy  content  for  both  upward 
and  downward  pulling  cases  at  <^o=70°  correspond  to  either  local  maxima  or  non- 
extrema.  Thus,  for  either  direction  of  puUing,  the  meniscus  shape  can  no  longer  be 
stable  for  such  high  aspect  ratios.  For  larger  h,/rb,  the  E  value  at  the  specified  0o 
is  a  local  maximum  of  the  energy  profile.  Hence,  all  such  menisci  are  not  stable  to 
small  disturbances,  and  tend  to  depart  from  the  given  meniscus  shape  and  contact 
angle.  Since  here  the  static  contact  angle  does  not  provide  a  minimum  energy  state, 
the  meniscus,  assuming  it  does  not  altogether  collapse,  will  tend  to  seek  such  a  state 
by  changing  shape  and  contact  angle  in  a  dynamic  manner, 

3.5.3 Ouasi-equilibrium  Motion  with  Hysteresis 

Contact  angle  hysteresis  is  generally  acknowledged  to  be  a  consequence  of 
three  factors:  (1)  surface  inhomogeneity,  (2)  surface  roughness,  and  (3)  impurities  on 
the  surface  (Miller  1985).  In  this  third  set  of  simulations  we  investigate  the  quasi- 


71 

equilibrium  path  negotiated  by  the  meniscus,  by  incorporating  a  simplified  hysteresis 
model  at  the  top  trijunction,  i.e., on  our  modeled  flat  plate.  Specifically,  the  top 
plate  is  oscillated  about  a  mean  height  h/  and  the  assumption  is  made  that  the  time 
scale  of  these  oscillations  is  such  that  a  succession  of  equilibrium  states  is  achieved. 
In  other  words,  the  Laplace- Young  equation  is  solved  for  each  position  h<.(t)  of  the 
oscillating  plate.  The  contact  line  at  the  die  is  allowed  to  remain  pinned  at  the  edge, 
as  long  as  the  Gibbs'  inequality  is  not  violated.  At  the  top,  the  boundary  condition 
is  modified  by  the  hysteresis  model  adopted  for  the  top  trijunction,  which  is  defined 
by  the  following  generally  observed  criteria: 

1.  The  contact  line  does  not  move,  i.e.,re(t)  is  a  constant  when  the  contact  angle 
varies  between  and  4>f^,  where  <i>^  and  are  respectively  the  advancing  and 
receding  contact  angles  (Dussan  V.  1979). 

2.  When  the  contact  angle  <^a  or  <^r  is  reached,  contact  line  motion  ensues.  To 
simplify  the  situation,  we  assume  that  the  contact  angle  does  not  vary  with  contact 
line  motion.  i.e.,<t>c=<i>^  or  <t>^  depending  on  the  direction  of  motion. 

3.  It  is  imposed  that  an  advancing  contact  line  cannot  recede  without  passing  through 
the  hysteresis  range  0^  to  0r.  Thus,  if  the  oscillatory  motion  of  the  plate  induces  an 
impending  recession  of  an  advancing  contact  line,  our  model  would  perforce  pin  the 
radius  r,(t)  at  such  a  location  and  require  0,  to  travel  from  0^  to  0r  before  allowing 
recession.   Similar  conditions  hold  for  the  receding  angle,  4>^. 

Schematic  presentation  of  the  hysteresis  model  is  shown  in  Figure  3-9.  An 
additional  simplification  resulting  in  the  model  is  relevant  in  the  presence  of  multiple 


72 

solutions.  When  the  contact  line  is  advancing,  for  instance,  we  choose  that  solution 
for  which  T^{t)  is  closest  in  the  direction  of  impending  motion.  So  also  for  a  receding 
contact  line.  Thus,  irrespective  of  the  stability  of  the  equilibrium  profiles,  proximity 
of  the  contact  line  is  made  the  criterion  for  choice  of  a  solution  from  those  available. 
Actually,  in  the  situation  where  dynamic  contact  angles  exist,  energy  criteria  in  terms 
of  the  existence  of  a  minimum  of  the  energy  curve  (  E  versus  )  is  no  longer 
applicable.  In  particular,  the  wetting  component  of  energy  E^^atd,  namely  TTj^{y,i-y,g) 
cannot  be  replaced  by  an  expression  such  as  7rre^7,gCos<^e.  for  otherwise  (7,i-7,g) 
would  cease  to  be  a  material  property.  Thus,  in  passing  through  the  successive 
equilibrium  states,  the  stability  of  menisci  cannot  be  deduced  in  terms  of  an  energy 
minimization  principle.  From  a  microscopic  viewpoint,  if  we  consider  that  the 
hysteresis  range  <t>j,  to  <^r  merely  represents  the  apparent  contact  angle  behavior,  and 
that  in  actuality,  the  microscopic  contact  angle  is  (f)„  (which  is  the  Young's  contact 
angle,  not  related  in  any  straightforward  way  to  0^  or  then  we  are  justified  in 
replacing  (7,r7sg)  by  7igCOS(^o-  However,  on  a  microscopic  scale,  the  phenomenon  of 
hysteresis  is  occasioned  by  the  presence  of  asperities  (Oliver  et  al.  1977,  Eick  et  al. 
1975)  and  surface  chemical  inhomogeneities  (Neumann  1972)  or  by  several  other 
effects  which  are  still  a  matter  of  investigation.  Thus,  the  wetted  energy  expression 
will  need  to  be  modified  to  7,gCos<j!)oA^  ,  where  A^^  will  incorporate  the  effects  of 
roughness,  etc.  No  tacit  assumption  regarding  E^,^  may  be  made  either  from  a 
microscopic  or  macroscopic  view,  and,  for  this  reason,  we  will  refrain  from  drawing 
any  conclusions  regarding  stability  of  the  quasi-equilibrium  path  simulated  hereunder. 
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It  may  be  pointed  out  that  the  situations  presented  in  this  section  and  the  two 
previous  ones  represent  two  limiting  types  of  behavior.  In  sections  3.5. land  3.5.2, 
the  boundary  conditions  and  stability  criterion,  in  terms  of  the  minimization  of  energy 
are  applicable  to  a  situation  where  the  time  scale  of  motion  of  the  top  contact  line 
is  greater  than  that  required  for  the  meniscus  to  achieve  equilibrium.  Here,  on  the 
other  hand,  it  is  supposed  that  the  time  scale  of  contact  line  motion  is  much  shorter 
than  that  required  for  reaching  an  equilibrium  state,  the  latter  being  determined  only 
by  incorporation  of  the  dynamics  of  the  interface  as  well  as  the  fluid  contained  in  the 
meniscus. 

Traditionally,  the  following  formula  has  been  employed  in  models  describing 
growth  of  a  crystal  via  a  pulling  process.  The  variation  of  radius  at  the  top  of  the 
meniscus  with  pull  rate  and  time  is  given  by  (Thomas  et  al.  1986,  Tatarchenko  and 
Brenner  1980), 

=       -  -^)  tan(4)  -  4)^  (3-13) 

where  (f>-<f>o  is  the  deviation  in  contact  angle  from  the  value  that  leads  to  fibre  growth 
with  non-uniform  diameter.  We  attempt,  in  what  follows,  to  evaluate  the 
applicability  of  this  expression  in  the  light  of  hysteresis.  In  the  event  of  0  being  a 
constant,  it  is  obvious  that  the  variation  of  the  radius  with  time  follows  the  same 
waveform  but  has  an  opposite  phase  with  respect  to  the  variation  of  meniscus  height, 
hc(t).  In  this  respect,  the  presence  of  hysteresis  introduces  a  non-linearity  in  the 
above  expression,  in  that  the  value  of  0  in  turn  depends  on  r<.(t)  and  its  rate  of 
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change.  However,  if  a  pinning  condition  arises,  i.e.,re  =  constant,  then  the  above 
formula  indicates  that,  provided  Up  ^  dh/dt,  <^  =  <Ao,  i.e.,  that  the  contact  angle 
remains  fixed.  But  this  is  neither  necessary  nor  true  when  hysteresis  is  included, 
because  it  can  happen  that  assumes  values  between  0a  and  0r,  while  the  contact 
line  remains  stationary.  The  results  of  our  calculations,  under  the  restriction  of 
equilibrium,  and  the  conditions  1^3,  prescribed  by  the  hysteresis  model  above,  will 
call  to  question  the  implications  of  Eq.  (3.13). 

Under  the  present  model,  the  oscillation  of  the  plate  is  effected  by  specifying 
h^(t)=he*+AhcSimrt,  where  in  all  the  results  presented,  Ahe=  0.2h/  was  given.  The 
hysteresis  range  0^  to  0r  is  also  specified,  with  the  values  corresponding  roughly  to 
the  material  properties  of  sapphire.  We  start  with  an  initial  condition  at  the  top 
injunction  which  facilitates  the  calculation  of  top  radius  re(t),  when  the  contact  angle 
is  (^^=<^^+A0A.  i.e-'the  contact  angle  is  assumed  to  have  been  displaced  so  that  the 
contact  line  motion  may  ensue.  But  for  the  initial  transients,  the  final  periodic 
behavior  of  the  meniscus  was  found  to  be  independent  of  the  specified  initial  state. 
The  subsequent  development  is  followed  by  plotting  the  time  variation  of  meniscus 
profiles,  top  radius  and  contact  angle.  As  in  the  previous  sections.  Bond  numbers  of 
10'^  and  2.5x10'^  are  used,  and  the  meniscus  is  internally  pressurized  to  a  given  ratio 
of  the  atmospheric  pressure.  Cases  simulated  below  relate  to  downward  pulling  of 
the  fibre. 

First,  a  low  Bond  number  case  (Bo  =10"*)  is  considered,  without  hysteresis,  the 
contact  angle  being  fixed  at  90°.  In  Figure  3-10,  the  sinusoidal  motion  of  the  plate 
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given  by  h,(t),  causes  a  corresponding  sinusoidal  motion  of  the  top  radius,  r<.(t).  The 
r,  variation  is  here  out  of  phase  with  h,(t)  as  represented  by  Eq.  (3.13).  It  is  noted 
that  the  radius  variation  in  Figure  3-10  for  the  imposed  amplitude  Ah„  is  quite  wide, 
assuming  values  from  0.62rbto  0.83rb.  Introduction  of  a  hysteresis  range  gives  rise, 
in  Figure  3-11,  to  a  square  waveform  for  r<.(t).  For  the  hysteresis  range  imposed, 
0R  =  85°  to  <f>f^=95°,  the  pinning  condition  comes  into  effect  when  4>^  lies  in  the 
hysteresis  range.  This  can  be  seen  from  Figs.  3-11  (c)  and  (d),  in  which,  while  <j)^ 
varies  between  4)^  and  0r,  the  radius  remains  a  constant.  It  is  not  difficult  to  foresee 
that  when  the  hysteresis  range  is  further  increased,  for  this  amplitude  of  oscillation, 
the  pinning  condition  is  eventually  enforced  throughout  the  cycle  of  hc(t),  while  (^..(t) 
assumes  a  sinusoidal  form  within  the  range  <^r  to  <^a-  In  terms  of  the  motion  of  the 
contact  line,  the  variation  of  r^H)  displayed  in  Figure  3-11  can  be  interpreted  as  a 
stick-slip  motion  (Jansons  1986). 

When  we  reach  a  higher  Bond  number,  say  2.5x10"^  as  in  Figure  3-12,  where 
the  contact  angle  is  fixed  at  100°  with  no  hysteresis  and  for  mean  height  h/  =  0.5, 
and  Ap  =5%.  Here  the  behavior  of  the  meniscus  registers  a  marked  change.  This 
is  consistent  with  the  difference  in  meniscus  behavior  arising  due  to  the  different 
Bond  number  regimes  as  seen  in  the  previous  sections.  It  is  particularly  noticed  here 
that  for  the  same  relative  variation  in  height  Ah,  as  in  Figure  3-10,  in  this  case  the 
variation  in  radius  r,  is  less  significant.  Also  the  radius  variation  is  in  fact  in  phase 
with  the  h,  variation,  being  significantly  in  opposition  to  Figure  3-10,  and  formula, 
Eq.  (3.13).  In  other  words,  the  radius  in  this  case  increases  with  height  instead  of 
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decreasing,  while  the  angle  remains  constant.  This  is  evidently  due  to  the  change  in 
meniscus  curvature  from  concave,  as  in  Figure  3-10,  to  convex.  For  a  higher 
meniscus,  h/  =  1.5,  for  the  same  4)^  and  Ap  shown  in  Figure  3-13,  the  radius 
variation  is  also  not  very  substantial  in  magnitude  and  the  waveform  is  not  identical 
to  h<.(t),  indicating  the  non-linearity  that  is  introduced  by  <^(t).  When  hysteresis  is 
present,  and  for  the  case  of  fairly  large  hysteresis,  and  tall  meniscus.  Figure  3-14 
shows  an  interesting  situation  arising  due  to  the  presence  of  multiple  solutions.  Note 
that  in  the  previous  cases,  only  one  solution  existed  in  each  case,  while  here  2  or  3 
solutions  are  obtained  for  certain  heights.  The  choice  of  solutions  is  made  as 
explained  above.  The  marked  excursions  of  radius  evident  in  the  figure  arise  when 
the  outer  solutions  seen  in  Figure  3- 14(a),  are  allowed  to  be  attained  by  the  imposed 
conditions  1^3.  The  hysteresis  model  in  this  case  causes  a  noticeable  change  in  the 
contact  angle  behavior  (Figure  3-14  (d))  in  comparison  to  the  static  model.  There 
appears  to  be  a  Bond  number  dependence  of  the  phase  relationship  between  drjdt 
and  dh,/dt,  which  is  not  reflected  in  the  relation,  Eq.  (3.13). 

Thus,  Eq.  (3.13)  relating  the  radius  with  height  of  the  meniscus,  seen  in  the 
light  of  a  quasi-equilibrium  path,  does  not  appear  to  exhibit  the  appropriate  features 
under  contact  angle  hysteresis.  In  the  absence  of  hysteresis,  the  Bond  number 
dependence  of  the  phase  is  not  captured.  When  the  Laplace- Young  equation  yields 
multiple  solutions,  the  choice  of  profiles  is  not  straightforward,  and  the  particular 
selection  procedure  adopted  indicates  an  abrupt  dynamic  behavior  which  may  not  be 
justifiably  represented  by  a  quasi-equilibrium  situation,  and  the  full  fluid-dynamic 
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consideration  is  required  to  assess  the  implication  of  such  behavior.  We  also 
reiterate  that  no  account  has  been  taken  of  the  stability  of  the  computed  equilibrium 
profiles. 

3.6  Conclusions 

1.  For  the  boundary  conditions  adopted  here,  which  pertain  some  crystal 
growth  processes,  namely  a  fixed  base  radius,  and  Neumann  boundary  condition  at 
the  top,  solutions  to  the  Laplace- Young  equation  were  obtained  using  a  fourth-order 
numerical  integration  scheme.  Multiple  solutions  were  obtained  for  a  given  contact 
angle  at  the  top. 

2.  The  energies  of  the  menisci  were  computed  and  stability  of  solutions 
decided  based  on  an  energy  minimization  concept.  While  in  many  cases  the  energy 
versus  4>^  curves  did  indeed  yield  a  minimum  at  a  specified  equilibrium  contact  angle 
00,  in  other  cases  local  maxima  or  non-extrema  were  obtained.  Depending  on  the 
operating  history  of  the  pulling  process,  it  may  transpire  that  an  unstable  meniscus 
may  form  first,  as  for  instance,  when  the  trijunction  point  of  an  unstable  shape  at  a 
new  hjTi,  is  closer  to  that  of  the  previous  equilibrium  position. 

3.  Unless  the  quantity  (h./rJ^Bo  is  large  enough,  hydrostatic  pressure  does  not 
considerably  influence  meniscus  shape,  and  the  effect  of  pulling  direction  is  therefore 
not  significant.  The  effect  of  external  pressurization,  Ap,  also  depends  on  the  Bond 
number.  For  higher.  Bond  numbers,  pressurization  can  markedly  widen  the  range  of 
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available  contact  angles.  Also,  for  sufficiently  tall  menisci,  pressure  (Ap)  and  gravity 
(Apgz)  interact  to  cause  menisci  to  assume  a  variety  of  curvatures. 

4.  For  modest  values  of  aspect  ratio  h,/rb,  there  always  exist  one  or  more 
stable  meniscus  shapes,  and  Young's  contact  condition  is  compatible  with  the 
requirement  of  energy  minimization.  For  higher  values  of  h^rb,  although  equilibrium 
solution  satisfying  Laplace- Young  equation  do  exist,  they  may  not  be  stable  solutions. 

5.  When  the  Young's  contact  condition  does  not  yield  a  minimum  in  energy 
content,  then  the  subsequent  non-equilibrium  dynamics  of  the  meniscus  needs  to  be 
simulated.  Such  treatments  have  been  reported  in  connection  with  other  situations 
(Young  and  Davis  1987,  Haley  and  Miskis  1981)  but  have  not  been  attempted  for  the 
pulling  process. 

6.  Quasi-equilibrium  studies  of  the  dynamics  of  menisci  have  been  made.  A 
simplified  hysteresis  model  was  employed  to  relate  the  contact  angle  at  the  top  to  the 
contact  line  motion.  The  conditions  enjoined  by  the  model  enforce  pinning  of  the 
contact  line  when  the  contact  angle  lies  between  0r  and  <f>^,  the  advancing  and 
receding  contact  angles.  It  is  found  that  the  relationship  in  terms  of  waveform  and 
phase  of  r,(t)  with  respect  to  h,(t)  cannot  be  captured  by  the  formula  usually 
employed  in  crystal  growth  simulations. 
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Figure  3-1  Schematic  of  meniscus  corresponding  to  edge-defmed  fibre  growth 
process. 


1 «  . 
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Figure  3-2     Characteristics  of  meniscus  formation  for  Bo  =  4x10"'  .  AP  =  0, 
hjT^  =0.5,  and  upward  puUingSO 

(a)  variation  of  <^),  with  angle  0^  . 

(b)  permissible  range  of  tJt^  with  respect  to  0,  . 

(c)  profile  shapes  for  different  angles  <i>^  . 

(d)  profiles  corresponding  to  0,,  =  60° . 

(e)  E  profile  v.s.0,  ,  for  the  specified  value  of  0„  =  60°  . 
(0  profiles  corresponding  to  <i>„  =  90° . 

(g)  E  profile  V.S.0,  ,  for  the  specified  value  of  0,  =  90°  . 
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Figure  3-2  --  continued 
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Figure  3-3  Characteristics  of  meniscus  formation  for  Bo=2.5xl0* ,  AP=5%, 
K^h  =2,  and  upward  pulling. 

(a)  profiles  corresponding  to      =  16°  . 

(b)  E  profile  v.s.d),  ,  for  the  specified  value  of  qi)„  =  16°  . 

(c)  profiles  corresponding  to  0„  =  20°  . 

(d)  E  profile  v.s.0,  ,  for  the  specified  value  of  0„  =  20°  . 
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Figure  3-4  Characteristics  of  meniscus  formation  for  Bo=2.5xI0-^  ,  AP=5% 
K^h  -2,  and  upward  pulling. 

(a)  profiles  corresponding  to  0^  =  60°  . 

(b)  E  profile  v.s.<t>^  ,  for  the  sp°ecified  value  of  0  =  60° 

(c)  profiles  corresponding  to  0^  =  90°  . 

(d)  E  profile  v.s.(^,  ,  for  the  specified  value  of  0  =  90° 
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Figure  3-5  Characteristics  of  meniscus  formation  for  Bo=4xlO"'  ,  AP=5%, 
hjTf,  =  5,  and  upward  pulling. 

(a)  profiles  corresponding  to      =  90°  . 

(b)  E  profile  v.s.0,  ,  for  the  specified  value  of  0^  =  90°  . 

(c)  profiles  corresponding  to      =  130°  . 

(d)  E  profile  v.s.0,  ,  for  the  specified  value  of  0„  =  130°  . 
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(i)  upward  pulling 


(ii)  downward  pulling 


Figure  3-6  Effect  of  the  direction  of  pulling  on  the  meniscus  characteristics 
with  Bo=10"^  and  AP  =  0. 

(a)  permissible  range  of  contact  angle. 

(b)  permissible  range  of  trijunction  location. 

(c)  permissible  meniscus  shapes  at  0„  =  10°  . 

(Note  that  with  given  Bond  number"  and  hJr,  less  than  1  the 
meniscus  shape  is  dominated  by  the  geometrical  constramt  from 
the  boundary  conditions.  .Hence  the  meniscus  shapes  are 
essentially  the  same  between  two  pulling  directions.  The  effect  of 
pulling  direction  becomes  noticeable  only  when  h  /r^  becomes 
larger.) 
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Figure  3-7  Effect  of  pressurization  on  the  meniscus  shapes  with  Bo  =  10"- and 
AP  =  1%.  (excess  interior  pressure). 

(a)  permissible  range  of  contact  angle. 

(b)  permissible  range  of  trijunction  location. 

(c)  permissible  meniscus  shapes  at      =  70°  . 

(Note  that  the  hydrostatic  pressure  effect  becomes  noticeable  only 
for  a  higher  h./r^  ,  as  evidenced  by  difference  caused  by  pulling 
directions.) 


87 


(a)  h^r,  =  0.5 


0.05 
0.045 

0.04 
0.035  i 

0.03  f 
E   0.025  \ 

0.02 
0.015 

0.01 
0.005 
0 


  opwtfd  puilins 

o  o  o  aownwara  pulling 

fii)  noiv<]nremuin 


I 


(i)  mininiism 


30  40  50  60  70  80  90  100110  120130 
<f>c 


(c)  hJT,  =  2.0 


40    45    50    53    60    65    70  75 


(b)  Vl,  =  1.0 


upwml  pulliOE 
>  downwara  pulimg 


(ii)  nofMxmmam 


(i)  minimain 
I 


20     30  40 


50  60 


70     80  90 


(d)  hJt,  =  2.5 


Figure  3-8  Energy  profiles  for  Bo  =  10"^  ,  AP  =  1%,       =  70°  and  various 
hJZf,  .  Both  upward  and  downward  pulling  cases  are  shown. 
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Meniscus  behavior  at  lower  Bond  number,  with  no  hvsteresis. 
r  Bo  =  10-^    AP  =  5%  h  ♦  =-  0  5  rA,^^_^jooo^.  ' 

(a)  Meniscus  profiles  obtained  for  the  duration  of  oscillation. 

(b)  Variation  of  height  h,  with  time. 

(c)  Computed  value  of  top  radius  r,  versus  time. 
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Fieure  3 


■11  Meniscus  behavior  at  lower  Bond  number,  with  small  hysteresis 

(  Bo  =  in-^    AP  ^        h' _^J]JUi„__SlZ^  =  85°) 

(a)  Meniscus  profiles  obtained  for  the  duration  "^niation 

(b)  Variation  of  height  h,  with  time. 

(c)  Computed  value  of  top  radius  r.  versus  time 

(d)  Computed  variation  of  contact  angle  0,  versus  time 
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Figure  3-12  Meniscus  behavior  at  higher  Bond  number,  with  no  hysteresis. 
(  Bo  =  2.5x10-'  ■  AP  =  5%.  h  '  =  0.5.6.  =       =  100°V 

(a)  Meniscus  profiles  obtained  for  the  duration  of  oscillation. 

(b)  Variation  of  height  h,  with  time. 

(c)  Computed  value  of  top  radius  r,  versus  time. 
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Meniscus  behavior  at  higher  Bond  number,  with  no  hysteresis 
(  Bo  =  2.5x10-   .AP  =  5%  h/  =  1.5.6.  =  6.  =  inn°S 

(a)  Memscus  profiles  obtained  for  the  duration  of  oscillation 

(b)  Variation  of  height  h,  with  time. 

(c)  Computed  value  of  top  Radius  r,  versus  time. 
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Figure  3-14  Memscus  behavior  at  higher  Bond  number,  with  hvsteresis 

f°  =  =  Jl.l^L5^,.^I10!^.^j00!}. 

a  Memscus  profiles  obtained  for  the  duration  of  oscillation 

(b)  Variation  of  height  h,  with  time. 

(c)  Computed  value  of  top  radius  r,  versus  time. 

(d)  Computed  variation  of  contact  angle  0,  versus  time 


Figure  3-15  Meniscus  behavior  at  higher  Bond  number,  with  small  hysteresis 
range. 

(  Bo  =  2.5x10-'  ■  AP  =  5%.  h,'  =  0.5.  6,  =  90°.  6^  =  88°'). 

(a)  Meniscus  profiles  obtained  for  the  duration  of  oscillation. 

(b)  Variation  of  height  h,  with  time. 

(c)  Computed  value  of  top  radius  r,  versus  time. 

(d)  Computed  variation  .of  contact  angle      versus  time. 


CHAPTER  4 

NUMERICAL  TECHNIQUES  FOR  SOLVING  FREE  AND  MOVING 
BOUNDARY  PROBLEMS  WITH  APPLICATION  TO  THE  EFG  PROCESS 


4.1  Introduction 

As  already  discussed,  prediction  of  the  interface  shape  and  movement  in 
the  EFG  process  is  possible  only  if  a  complete  dynamic  model  incorporating  all 
the  geometric  and  thermal  aspects  of  the  growth  process  is  available.  In  order  to 
facilitate  such  simulations,  various  computational  elements  must  be  developed, 
including  appropriate  time  stepping  and  spatial  discretization  schemes,  thermal 
field  solutions  in  the  irregularly  shaped  domain  by  the  melt/gas  meniscus  and  the 
melt/solid  interface  and  method  for  tracking  the  free  and  moving  boundaries 
(Tseng  et  al.  1992,  Imaishi  et  al.  1991,  VoUer  et  al.  1990,  Floryan  and  Rasmussen 
1989,  Brown  1988,  Salcudean  and  Abdullah  1988,  Viskanta  1988,  1985  and  1983, 
Crank  1984,  Ungar  and  Brown  1984,  Epstein  1983). 

A  literature  survey  of  some  existing  methods  for  solving  the  Stefan 
problems  is  given  in  section  4.2.  Then  the  issues  arising  from  developing  the 
front  tracking  techniques  are  investigated  in  section  4.3.  We  assess  the 
applicability  of  the  oft-used  quasi-stationary  approximation  for  both  small  Stefan 
number  and  large  Stefan  number  cases.  A  new  generic  interface  tracking 
procedure,  based  on  a  combined  L^grangian/Eulerian    framework,  is  designed; 
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this  procedure  has  been  developed  in  conjunction  with  the  previously  mentioned 
time  stepping  scheme  and  the  meniscus  model  to  simulate  the  EFG  process. 

After  the  presentation  of  the  necessary  numerical  techniques,  the  EFG 
process  physics  is  described  in  section  4.4.  Fluid,  thermal,  and  geometrical 
related  non-dimensional  parameters  in  thin  crystal  growth  via  the  EFG 
technique  are  evaluated.   Based  on  the  order  of  magnitude  estimation,  a 
mathematical  model  is  developed  for  the  energy  transport  of  the  EFG  process  in 
section  4.5.  Combining  the  front  tracking  and  moving  gridding  techniques,  time 
stepping  scheme,  solution  procedure  of  the  meniscus  shape  and  contact 
condition,  and  the  detailed  thermal  fields  in  both  liquid  and  solid  phases,  a 
complete  numerical  tool  for  the  EFG  process  is  presented  in  section  4.6.  The 
model  is  then  applied  to  analyze  cases  with  two  Stefan  numbers,  i.e., St  =  1  and 
0.024.  The  steady  state  solutions  corresponding  to  the  specified  boundary 
condition  are  obtained  for  both  cases.  It  is  identified  that  an  appropriate 
scaling  which  accounts  for  the  relative  magnitude  of  the  latent  heat  is  necessary 
to  properly  capture  the  interface  movement  and  to  improve  the  computational 
efficiency.  Based  on  the  computed  results  of  the  two  cases  with  different  values 
of  Stefan  number,  some  conclusions  for  steady  state  solutions  are  made  in 
section  4.7.  The  obtained  steady-state  solutions  will  be  used  as  base  solutions  to 
study  the  dynamic  behaviors  of  EFG  system  in  Chapter  5. 
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4.2  Literature  Review 
The  characteristic  feature  of  any  phase  change  problem  is  the  coupling  of 
thermal  fields  in  the  different  phases  with  free  and  moving  boundaries  that  not  only 
separate  each  phase,  but  also  dynamically  evolve.  The  free  surface  and  propagating 
phase  front  make  such  problems  nonlinear.   Because  of  the  nonlinearity,  only  a  few 
exact  solutions  have  been  developed,  and  these  are  limited  to  simple  geometries 
and  boundary  and  initial  conditions.  For  most  practical  cases,  the  solutions  are 
obtained  using  numerical  methods.   In  what  follows  a  brief  survey  of  the  existing 
methods  for  solving  the  phase  change  problems  is  presented. 

4.2. 1  Exact  Solutions 

One-dimensional  heat  transfer  with  phase  change  was  first  posed 
mathematically  by  Stefan  and  later  solved  analytically  by  Neumann  (Crank  1984, 
Ozisik  1980,  Carslaw  and  Jaeger  1959).  The  exact  solutions  are  limited  to  a 
number  of  idealized  situations  involving  semi-infmite  regions  or  infmite  regions 
subject  to  simple  boundary  and  initial  conditions.  To  date,  no  exact  solutions  have 
been  developed  for  phase  change  problems  where  the  phase  change  front  is  a 
function  of  multiple  space  coordinates. 

4.2.2  Approximate  Solutions 

The  integral  method,  which  dates  back  to  von  Karman  and  Pohlhausen, 
who  used  it  for  the  approximate  analysis  of  boundary  layer  equations,  was  applied 
by  Goodman  (1958)  to  solve  a  one-dimensional  transient  phase-change  problems. 
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This  method  provides  a  relatively  straightforward  and  simple  approach  for 
approximate  analysis  of  one-dimensional  transient  phase-change  problems. 

4.2.3  Numerical  Solutions 

Many  kinds  of  crystals  grown  by  a  EFG  technique  usually  have  high 
melting  temperatures;  for  example,  sapphire  (AI2O3)  with  a  melting  temperature 
of  1046  K.  Therefore,  the  need  of  simulating  the  solidification  process  has  led  to 
a  development  of  several  numerical  methods.  Most  of  them  were  summarized  in 
Crank  (1984),  Wilson  et  al.  (1978),  and  Ockendon  and  Hodgkins  (1975). 

Chemouskou  (1970)  proposed  the  isotherm  migration  method,  which 
interchanges  the  dependent  variable  with  the  space  variable  for  one-dimensional 
phase  change  problems.  Crank  and  Gupta  (1975)  were  the  first  to  apply  the 
isotherm  migration  method  to  phase  change  problems  in  two  dimensions.  The 
advantages  of  this  technique,  in  tracking  a  moving  isotherm  at  the  fusion 
temperature,  are  obvious,  since  the  phase  boundary  is  itself  an  isotherm,  provided 
the  phase  change  takes  place  at  a  constant  temperature.    Results  predicted  by  this 
approach  have  been  found  to  be  satisfactory.  Saitoh  (1978)  extended  the  two- 
dimensional  isotherm  migration  method  from  regular  two-dimensional  geometries 
to  arbitrary  shaped,  doubly  connected  two-dimensional  geometries. 

The  moving  heat  source  (or  the  integral  equation)  method,  originally 
applied  by  Lightfoot  (1929),  is  based  on  the  concept  of  representing  the  liberation 
(or  absorption)  of  latent  heat  by  a  moving  plane  heat  source  (or  sink)  located  at 
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the  solid/liquid  interface.  The  temperature  at  any  point  will  be  due  to  the 
superposed  contributions  from  this  moving  source  (or  sink).  In  this  method,  the 
analysis  of  the  phase-change  problem  is  transformed  to  the  solution  of  an  integral 
equation  for  the  location  of  the  liquid/solid  interface.   Ozisik  (1978)  discussed  the 
use  of  a  moving  heat  source  in  the  conduction  equations  to  account  for  the  latent 
heat  effect.  He  presented  a  mathematical  development  that  explicitly  casts  the 
moving  boundary  problem  into  a  standard  heat  conduction  problem  with  a  moving 
heat  source.  By  doing  this,  the  solution  can  immediately  be  written  in  terms  of 
Green's  function,  and  numerically  implemented. 

In  general  the  numerical  methods  for  solving  the  phase  change  problems 
can  be  divided  into  two  categories:  the  fixed-grid  approach  and  the  transformed- 
grid  approach.   The  fixed-grid  approach  implies  that  a  grid  is  fixed  in  space  and 
the  interface  conditions  are  accounted  for  by  the  definition  of  suitable  source 
terms  in  the  governing  equations.   In  the  transformed-grid  approach,  the 
governing  equations  and  their  boundary  conditions  are  cast  into  a  generalized 
curvilinear  coordinate  system,  so  the  grid  might  adapt  with  the  moving  freezing 
front.  A  comparison  between  transformed  grids  and  fixed  grids  on  a  test  problem 
of  the  melting  of  gallium  was  made  by  Lacroix  and  VoUer  (1990).  They  have 
found  that  if  an  efficient  grid  generation  can  be  achieved,  the  CPU  usage  of  the 
transformed  grid  is  very  close  to  that  of  the  fixed  grid.  Furthermore,  a  fixed  grid 
can  produce  accurate  predictions  with  the  same  order  of  mesh  size  as  that  used  by 
a  transformed  grid.  Thus,  according  to  these  authors,  no  conclusive  preference 
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can  be  made  for  either  method.   However,  the  tests  covered  only  a  limited  range 
of  grid  and  time  resolution;  more  effort  in  this  direction  will  be  needed  to  further 
clarify  this  issue. 

One  well  known  method  to  deal  with  phase  change  problems  on  fixed  grids 
is  the  enthalpy-porosity  technique.   This  approach  easily  permits  solutions  for 
substances  whose  change  of  phase  occurs  over  a  range  of  temperature.    Brent  et 
al.  (1988)  have  adopted  this  technique  in  conjunction  with  the  SIMPLE  procedure 
proposed  by  Patankar  (1980)  to  simulate  the  two-dimensional  convective  melting 
of  pure  gallium  in  a  rectangular  cavity.  The  enthalpy-porosity  technique  was 
validated  by  comparing  its  results  with  the  experimental  results  obtained  by  Gau 
and  Viskanta  (1986).  A  review  of  available  implicit  difference  enthalpy-porosity 
schemes  can  be  found  in  Voller  (1985). 

One  of  the  drawbacks  of  the  enthalpy  formulation,  particularly  for  the 
case  where  the  change  of  phase  occurs  at  a  single  temperature,  is  that  a  plot  of 
temperature  against  time  for  a  given  node  tends  to  exhibit  a  "plateauing" 
tendency.  This  arises  out  of  the  fact  that  a  computational  grid  models  or 
represents  a  discrete  region  in  space.  Obviously,  it  requires  a  finite  amount  of 
time  to  melt  a  discrete  region.  As  a  consequence  of  a  node  being  held  fixed  at  a 
single  temperature  for  a  discrete  amount  of  time,  the  effect  is  also  felt  by  the 
surrounding  nodes,  causing  a  plateauing  of  the  temperature.    Voller  and  Cross 
(1981)  have  developed  a  smoothing  technique  that  can  be  applied  to  a  final  set  of 
numerical  results.  This  smoothing  technique  has  the  effect  of  bringing  the  time- 
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temperature  history  of  a  given  node  into  better  agreement  with  other  published 
results. 

Treatments  of  the  phase  change  and  the  associated  transport  process  based 
on  the  fixed  grid  method  have  been  developed  to  treat  quite  complicated 
problems.  Bennon  and  Incropera  (1987)  have  used  a  volume  averaged  continuum 
model  with  the  SIMPLE  algorithm  (Patankar  1980)  to  investigate  solidification  of 
binary,  aqueous  ammonium  chloride  solution  in  a  rectangular  cavity.  Beckermann 
and  Viskanta  (1988)  have  also  used  the  same  numerical  method  to  solve  the 
complex  phase  change  problem.  Shyyand  Chen  (1991,1990)  have  combined  the 
enthalpy  formulation  and  computational  scheme  based  on  an  adaptive  grid  in 
curvilinear  coordinates  to  solve  the  phase  change  problems  with  high  Rayleigh 
and  Marangoni  numbers,  and  under  both  normal  and  reduced  gravity  conditions. 
With  relatively  finer  and  improved  resolutions,  various  detailed  transport 
phenomena  and  interface  shapes  have  been  reported.   Recently,  a  turbulence 
model  has  been  incorporated  into  this  framework  to  predict  the  solidification 
aspects  of  a  continuous  ingot  casting  process  (Shyy  et  al.  1992a),  yielding  a  highly 
advanced  tool  for  analyzing  transport  processes  during  solidification. 

Crowley  (1983)  applied  an  explicit  finite  difference  scheme  to  an  enthalpy 
formulation  and  coordinate  transformation  techniques  to  study  the  crystal  pulling 
by  the  Czochralski  technique.   She  applied  her  model  to  analyze  the  initially 
steady  Czochralski  system  with  a  base  temperature  perturbed  at  the  melt  inlet. 
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Based  on  the  computed  response  of  trijunction  location,  she  concluded  that  the 
simulated  configuration  is  unstable  to  such  disturbances. 

A  numerical  solution  for  melting  around  a  vertically  heated  cylinder 
incorporating  the  effects  of  natural  convection  was  obtained  by  Sparrow  et  al. 
(1977).  They  used  the  coordinate  transformation  to  immobilize  the  solid-liquid 
interface  and,  consequently,  appear  to  have  been  the  first  investigators  to  apply  it 
to  a  multidimensional  convection-dominated   Stefan  problem.   However,  the 
higher  order  curvature  terms  for  the  interface  were  neglected,  limiting  the  validity 
of  the  method  to  cases  where  the  radius  of  the  interface  varies  only  slightly. 
Furthermore,  a  quasi-steady  assumption  was  invoked  where  the  effects  of 
interface  motion  were  ignored;  hence  the  method  can  treat  the  problem  only  with 
very  modest  phase  change  rate.  Thompson  and  Szekely  (1989)  have  applied  this 
quasi-steady  approach  to  predict  the  phase  change  problem  in  a  cavity. 

The  above  mentioned  works  are  based  on  the  finite  volume/difference 
formulations.   Silliman  and  Scriven  (1980)  and  Cuvelier  and  Schulkes  (1990)  have 
used  the  finite  element  method  (FEM)  to  calculate  the  predictions  of  steady 
viscous  free  surface  flows.  In  their  analysis,  only  liquid/gas  interface  is 
considered,  and  phase  change  is  not  included.  Albert  and  O'Neill  (1986)  and 
Lynch  (1982)  used  a  finite  element  method  with  deforming  elements  to  track  a 
moving  phase  front  and  to  model  phase  change  problems. 

Ettouney  and  Brown  (1983)  have  used  a  finite  element  analysis  to  model 
the  steady  EFG  process.  Several  variations  of  the  methods  have  been  developed 
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by  this  group  (Sackinger  et  al.  1989,  Derby  and  Brown  1986),  depending  on  the 
choice  of  coordinates  (original  or  transformed),  distinguished  boundary  condition 
(local  heat  balance  or  thermal  equilibrium),  and  the  method  of  solving  the 
resulting  system  of  equations  (successive  approximation  method  or  Newton's 
method).   Their  numerical  method  was  limited  to  the  steady  solidification 
problems. 

Crochet  et  al.  (1987),  Dupont  et  al.  (1987)  and  Wouters  et  al.  (1987) 
investigated  the  oscillatory  fluid  flow  in  a  horizontal  Bridgman  technique.  The 
method  was  based  on  an  implicit  time-marching  technique  using  FEM  but  was  not 
applied  to  the  moving  solid-liquid  interface  problem.   Other  works  in  this  area 
can  be  found  in  (Heinrich  et  al.  1989,  and  Dantzig  1989). 

In  order  to  accurately  resolve  the  solid/liquid  interface  for  the  isothermal 
phase  change  problem  encountered  here,  the  explicit  front  tracking  method  is 
adopted.   A  new  technique  has  been  developed,  as  presented  below. 

4.3  Front  Tracking  and  Moving  Grid 
The  focus  of  this  section  is  the  numerical  simulation  of  interface  motion 
during  melting/solidification  of  pure  material.   First,  we  assess  the  applicability  of 
the  aforementioned  quasi-stationary  approximation  (Lacroix  1989,  Thompson  and 
Szekely  1989,  Crowley  1983,  Sparrow  et  al.  1977)  for  tracking  the  interface 
motion.  Such  an  approximation  for  small  St  as  well  as  large  St  is  numerically 
examined.  Next,  a  generic  interface  tracking  procedure  is  designed,  which 
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overcomes  restrictions  of  single-validness  of  the  interface  imposed  by  commonly 
used  mapping  methods  (Lacroix  1989,  Thompson  and  Szekely  1989,  Crowley  1983, 
Sparrow  et  al.  1977).  This  method  incorporates  with  ease  interface  phenomena 
involving  curvature,  which  assumes  importance  at  the  smaller  scales  of  a  deformed 
interface. 

Among  the  methods  proposed  in  the  literature  to  track  moving  fronts 
(Floryan  and  Rasmussen  1989,  Crank  1984)  one  commonly  used  approach 
(Lacroix  1989,  Thompson  and  Szekely  1989,  Crowley  1983,  Sparrow  et  al.  1977) 
involves  coordinate  transformation  to  map  the  irregularly  shaped  physical  domain 
onto  a  regular  computational  domain.  The  solutions  may  then  be  obtained  in  the 
transformed  domain  and  the  interface  position  updated.  Often  a  further 
simplification  is  made  by  assuming  that  the  time  scale  for  interfacial  motion  is 
long  compared  to  thermal  relaxation  times.  This  justifies  dropping  terms 
representing  grid  movement  from  the  equations  governing  interface  motion.  This 
procedure  is  termed  quasi- stationary  and  is  appropriate  for  slow-moving 
interfaces. 

4.3.1  Assessment  of  the  Quasi-stationary  Approximation 

As  a  result  of  the  nonlinearity  of  the  equation  introduced  by  the  motion  of 
the  phase  change  boundary  there  are  but  a  few  exact  solutions  of  problems 
concerned  with  melting  and  solidification  (Crank  1984,  Ozisik  1980,  Carslaw  and 
Jaeger  1959).  One  such  problem  is  the  melting  of  a  solid  (Fig.  4-1)  confined  to  a 
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semi-infinite  region.  The  solution  is  known  as  Neumann's  solution  (Crank  1984. 
Ozisik  1980,  Carslaw  and  Jaeger  1959).  We  will  use  it  to  verify  the  developed 
numerical  algorithm. 

For  a  pure  conduction  problem  with  phase  change  between  liquid  and  solid, 
the  governing  equations,  in  dimensional  form,  for  the  energy  transfer  and  interface 
movement  are  respectively  (Crank  1984) 


The  partials  with  respect  to  n  represent  derivatives  in  the  direction  of  the  local 
normal  to  the  interface.  The  following  non-dimensionalizing  procedure  is  adopted. 
X=x/lr  ,    Y=y/lr  ,  S=s/lr  ,  N=n/lr  ,  1^  being  a  representative  length  scale  to  be 
chosen  for  the  specific  problem  studied.  r'=  ol-^\^,    p*  =  Pj  /pi  and  6,  =  (Tj  - 
Tm)/(T,  -  T^  ,  where  Tn,  and  T,  are,  respectively,  the  melting  temperature  and  the 
imposed  temperature  at  the  appropriate  boundary.  Consider  heat  transport  in  the 
liquid  phase  only,  i.e.,^,  =  0.    Then,  Eqs.  (4.1)  and  (4.2),  discarding  the 
superscript,  become: 


i  =  l,s 


(4.1) 


(4.2) 
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dd    _  d^e    .    d\  (4j) 


+ 


dr'      dX^  BY^ 
_  dd  _  p  *  dS  (4.4) 

'dN  IFTr 

where  St  =  Cp(  T„-TJ  /If  is  the  Stefan  number. 

If  X=S(Y,t),  then  Eq.  (4.4)  specialized  to  an  isothermal  interface  can  be 
rewritten  as 

—  =  -  £l—  ^  (4.5) 

dY      dX  St  dr' 

Following  the  transformation  from  Cartesian(X,Y,r)  to  curvilinear  (^,v,t') 

coordinates  (Thompson  et  al.  1985),  such  that  X  =  X  (|,r/,T*),Y  =  Y{^,7],t)  and 

r*  =  T,  Eqs.  (4.3)  and  (4.4)  become: 


d  +  1(X  Y  -YX  )d,  +  -(Y^  -XJ  )d 


(4.6  a) 


^[(Ye\  -  (Y^d)^  [1  ^  (1[-(X5),  ^  iX^S)^y  ] 

=  -A[^+I(xy  -  YX)s,  +  l(y^  -  A:,y)5j 
St  dr  J  "  '      "         y  ^ '      ^  ' 


(4.6  b) 


where 


(4.7) 
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Eqs.  (4.6  a)  and  (4.6  b)  constitute  the  complete  set  of  the  equations  in  the 


generalized  coordinates  for  energy  transfer  and  interface  movement.   In  the 
course  of  the  computations,  at  every  time  step,  both  Eq.(4.6  a)  and  (4.6  b)  need 
to  be  iteratively  updated  until  they  are  simultaneously  satisfied  everywhere  in  the 
domain.  Standard  second-order  central  differences  are  used  for  all  spatial 
derivatives.  The  quasi-stationary  approach  adopted  by  many  researchers,  e.g., 
Lacroix  (1989),  Thompson  and  Szekely  (1989),  Crowley  (1983),  Sparrow  et  al. 
(1977),  simplifies  the  above  equations  by  dropping  all  the  terms  involving 
coordinate  movement,  i.e.,X„  Y^.  To  test  whether  such  a  simplification  is 
acceptable,  we  design  a  situation  where  the  interface  moves  as  a  planar  front.  In 
addition,  the  temperature  in  the  solid  is  considered  uniform,  i.e., 6^  =  0.  Thus  the 
heat  transport  is  effectively  one-dimensional  and  occurs  in  the  liquid  phase  only. 
Boundary  and  initial  conditions  are  therefore  given  by  Crank  (1984), 


B.C. 


e 


1,  X=0,  T>0 


(4.8) 


e 


0,  X>5,  r>0 


I.e. 


erfi 


X 


(4.9) 


5(X,r„=0.1)  =  2\f^, 
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where  X  is  the  root  of  the  following  equation 

Xe'^'em  =  A  (4.10) 

The  exact  solution  (Crank  1984)  for  Eqs.(4.3)  and  (4.4),  with  p*  =  1,  with 
the  boundary  and  initial  conditions  given  in  Eqs.(4.8)  and  (4.9)  is 

'^77=^  (4.11) 

em 

Two  numerical  solution  procedures  are  designed  to  compare  the  solutions 
obtained  from  the  complete  form  of  Eqs.  (4.6  a)  and  (4.6  b)  and  the  quasi- 
stationary  form.  These  are, 

(1)  Full  treatment:  A  standard  procedure  involving  backward  Euler  time  stepping 
along  with  a  second-order  central  difference  spatial  discretization  scheme  is  used 
to  solve  Eqs.  (4.6  a)  and  (4.6  b).  This  full  set  of  equations  is  solved  iteratively  in  a 
coupled  manner  to  continually  update  the  nonlinear  coefficients  resulting  from  the 
coordinate  movement  and  transformation. 

(2)  Quasi-stationary  treatment:  By  invoking  the  quasi-stationary  assumption, 
coordinate  movement  terms  are  neglected.  On  dropping  the  grid  movement  terms, 
the  equations  governing  energy  transfer  and  interface  movement  are  decoupled 
and  only  a  simple  explicit  procedure  is  needed  to  update  the  interface  location. 

Three  values  of  the  Stefan  number  were  chosen  to  test  the  performance  of 
the  above  two  numerical  approaches.    Eleven  grid  points  are  used  along  the  x- 
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direction  in  each  case.  For  small  St.  i.e.,0.1303,as  seen  in  Fig.  4-2  (a),  both  full 
and  quasi-stationary  solutions  are  in  close  agreement  with  the  exact  solution  in 
terms  of  the  interface  trajectory  as  well  as  the  temperature  profile,  although  the 
full  approach  is  marginally  superior.  However,  as  expected,  with  increasing  St,  the 
quality  of  solutions  obtained  from  the  quasi-stationary  method  is  progressively 
degraded.   In  contrast,  good  agreement  is  maintained  between  the  solutions  of  the 
full  approach  and  the  exact  solution.  Obviously,  with  a  low  St,  the  interface 
velocity  is  modest,  and  neglecting  grid  movement  terms  does  not  impact 
significantiy  on  the  numerical  accuracy.  As  St  becomes  larger,  however,  the 
simplification  to  the  quasi-stationary  approach  is  no  longer  acceptable. 

4.3.2  A  General  Procedure  for  Interface  Tracicing 

As  mentioned  in  previous  subsection,  Eqs.  (4.6  a)  and  (4.6  b)  apply  only  to 
an  isothermal  interface.   Also,  the  construction  of  those  equations  precludes  the 
possibility  of  capturing  branched  interfaces.  As  a  first  step  towards  overcoming 
this  difficulty,  an  alternate  interface  treatment  is  designed  as  follows.  Consider 
the  equation  for  interface  advance 


k  dd 

s  s 


k,  BN 


(4.12) 
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The  local  normal  to  the  interface  is  given  by 


(4.13) 


N  = 


|VF|  dx 


where  F  =  F(X,Y,t')  is  the  curve  defining  the  interface.  The  interface  shape  is 
defined  in  a  piecewise  fashion  to  facilitate  handling  of  branched  interfaces.  Here 
a  quadratic  polynomial  fit  is  performed  for  three  successive  nodal  points  at  each 
point  of  the  interface.  Thus,  at  the  i*  point  on  the  interface  we  designate  the 
curve,  Yi  =  a^X;^  +  b^X^  +  q  ,  i.e.,Fi  =  Yi  -  (  a^Xj'  +  bjX;  +  C; )  defines  the 
interface.  The  a^,  bj  and  Cj  are  determined  from  the  known  values  of  (Xj,Yj),  j  = 
i-l,i  ,  i+1.  The  choice  of  the  piecewise  approximating  polynomial  is  not 
restricted  to  the  parabolic  form.  In  fact,  using  circular  arc  elements  is  more 
convenient  for  handling  multiple-valued  interfaces  and  is  in  use  in  ongoing  work. 
The  local  curve  definition  yields  the  derivatives  (F,,  Fy  and  F^  )  at  each  point  on 
the  interface.   We  may  write  Eq.(4.12)  as 


In  computing  the  interface  normal  velocity  then,  one  seeks  to  obtain  the 
derivatives  F^,  Fy  and  6^  ,  dy.  The  derivatives  of  temperature  may  be  obtained  in 
the  transformed  coordinates  itself.  F^  and  Fy  of  course  are  directly  available  in 
the  physical  domain  from  the  curve  fit.  Thus  the  new  coordinates  of  the 
interfacial  points  are  obtained  from: 


|VF| 


[  (FA  -  F^^),  -  (FA  -  F^,).  ] 


(4.14) 


Ill 

Xn^i  ^X"  ^  —-!^8t  (4.15  a) 

dX  |VF| 

y*'  =  y-  +  K-!^8t  (4.15  b) 

BY  |VF| 

where  5r  is  the  non-dimensional  time  step  size.  Having  obtained  these  new 
coordinates  of  the  curve,  the  thermal  field  is  solved  for  once  again,  the  curve  fit  is 

-  h  .... 

performed  at  the  interface  and  a  fresh  interface  position  is  obtained  from  Eqs.(4.15 
a)  and  (4.15  b).  All  these  procedures  are  performed  in  a  fully  coupled  manner 
involving  interaction  among  the  temperature  field,  interface  motion  and  grid 
movement  at  each  iteration.  This  alternative  interface  ti^eatment  is  compared  first 
with  the  results  of  section  4.2.1, where  the  thermal  field  is  active  in  the  liquid 
region  only.  The  boundary  and  initial  conditions  are  given  by  Eqs.  (4.8)  and  (4.9). 
From  Fig.  4-3,  it  is  evident  that  for  all  three  Stefan  numbers  tested,  namely,  0.1303, 
1.2216  and  2.8576,  the  accuracy  of  the  present  scheme  is  satisfactory  and 
comparable  to  the  full  treatment  in  section  4.2.1, Fig.  4-2.  As  can  be  seen,  the 
computed  and  exact  temperature  fields  are  in  excellent  agreement. 

The  metiiod  has  been  extended  to  the  case  where  temperature  gradients 
exist  in  both  liquid  and  solid,  and  applied  to  study  the  development  of  a 
morphologically  unstable  interface  (Shyy  et  al.  1992  b). 

4.3.3  Assessment  of  Time-stepping  Schemes  for  Phase  Change  Problems 

As  already  mentioned,  the  EFG  process  operates  in  the  conduction 
dominated  regime.  Hence,  the  interaction  between  the  time  stepping  and  the 
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convection  schemes  discussed  in  Chapter  2,  while  generically  important,  is  not 
critical  for  our  main  interest  here.  In  what  follows,  a  comparative  study  is 
conducted  to  assess  the  relative  performance  between  the  backward  Euler  and 
Crank-Nicolson  schemes  for  solving  one-dimensional  Stefan  problems. 

The  governing  equations  of  the  unsteady  heat  conduction  problem  in  one 
phase,  along  with  the  energy  balance  requirement  at  the  interface  boundary, 
boundary  conditions  and  initial  condition  are  given  by  Eqs.  (4.3),  (4.4),  (4.8)  and 
(4.9). 

The  computations  start  with  the  exact  solution  at  t  =  0. 1  as  the  initial 
condition,  and  stop  at  r  =  0.2  with  four  different  values  of  At  with  St  =  0.1303. 
The  predicted  interface  locations  at  r  =  0.2  along  with  required  CPU  time  using 
Micro Vax  3100  are  listed  in  Tables  4-1  (a)  and  4-1  (b).  Two  different  spatial  grid 
sizes.  Ax  =  0.0158  with  11  grid  points  and  Ax  =  0.0079  with  21  grid  points,  along 
with  several  time  step  sizes  have  been  used  in  the  assessment. 

Based  on  the  Table  4-1  (a)  the  accuracy  of  the  Crank-Nicolson  method  is 
very  comparable  to  the  backward  Euler  method.   In  terms  of  the  CPU  time 
requirement,  both  methods  are  essentially  the  same,  as  expected.  It  should  be 
pointed  out  that  in  the  interface  updating  procedure  via  the  integration  of  the 
paths  of  marker  particles,  discussed  in  section  4. 3. 2, a  first-order  one-sided  scheme 
is  used  to  estimate  the  thermal  gradients  at  the  interface.   Consequently,  the 
marker  velocity  is  first-order  accurate,  and  the  resulting  interface  location  is  not 
dictated  by  the  choice  of  time  stepping  scheme.  It  is  apparently  this  practice  that 
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is  responsible  for  the  comparable  performance  between  backward  Euler  and 
Crank-Nicolson,   However,  as  will  become  clear,  due  to  the  nonlinearity  of  the 
phase  change  problem,  only  a  relatively  small  time  step  size,  around  10'\  can  be 
used  to  sustain  a  stable  computation.    Accordingly,  there  is  little  preference 
between  the  two  choices  of  time  stepping  schemes.  For  all  the  computations  to 
be  presented  in  the  following,  the  backward  Euler  scheme  has  been  adopted. 

4.4  EFG  Process  Physics 

In  the  following,  the  basic  process  physics  will  be  first  introduced  along 
with  assessments  of  the  key  control  parameters  including  the  Fourier,  Bond, 
Rayleigh,  Marangoni,  Peclet  numbers  and  aspect  ratio.  Then  the  appropriate 
governing  equations  and  boundary  conditions  will  be  presented,  and  the  numerical 
aspects  needed  for  simulating  such  a  process  discussed  in  the  next  section. 

The  basic  EFG  operation  considered  here  is  shown  in  Figures  1-1  and  1-2, 
including  a  schematic  illustration  and  a  description  of  the  boundary  conditions 
and  notation  adopted.   A  major  feature  of  the  present  process  is  the  relatively 
small  size  of  the  meniscus  and  the  fibre  diameter.   As  indicated  by  the  Bond 
number,  defined  as 

Bo  =  (4.16) 

where  rb  is  the  lower  contact  radius  on  the  die,  7,^  is  the  surface  tension  between 
the  melt  and  the  air,  Ap  is  the  density  difference  between  the  melt  and  the  air. 
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and  g  is  the  gravitational  acceleration.    With  the  size  and  properties  given  in 
Table  4-2,  Bo  is  of  the  order  of  10"*,  meaning  that  the  capillary  effect  is  important 
in  determining  the  meniscus  shape.  Three  types  of  convection  exist  in  the  melt  in 
such  a  system.  The  first  one  is  the  buoyancy  driven  convection  caused  by  the 
interaction  between  density  gradient  and  gravitational  acceleration  in  the  bulk  of 
the  melt.  The  second  one  is  the  thermocapillary  convection  caused  by  the  surface 
tension  gradient  due  to  the  temperature  gradient  along  the  melt/gas  interface. 
The  third  one  is  forced  convection  caused  by  the  flow  of  the  melt  towards  the 
solidification  front,  due  to  the  crystal  pulling,  in  order  to  conserve  mass.  Based 
on  the  nominal  values  fibre  diameter  of  0.005  ~  0.025  cm,  and  a  temperature 
variation  of  5  ~  25  K  from  the  die  tip  to  the  solid-liquid  interface,  the  Rayleigh 
number 


is  0(10'*)  for  sapphire,  according  to  the  property  values  listed  in  Table  4-2.  Since 
the  interface  separating  the  melt  from  the  outside  gaseous  environment  involves  a 
large  temperature  gradient  along  the  free  boundary,  the  Marangoni  effect  should 
also  be  estimated.    Again,  with  the  reference  values  given  in  Table  4-2  and  a 
surface  tension  of  -  0.06  dyne/cm-K  for  dy^/dT,  the  Marangoni  number 


Ra  = 


(4.17) 


Ma  = 


(4.18) 


is  about  0(1).  The  Peclet  number  measures  the  importance  of  heat  transfer  by 
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convection,  based  on  the  pull  speed,  Up,  in  the  growth  direction  relative  to 
conduction 


Pe  =  ^  (4.19) 

is  0(10"').  For  the  present  system,  the  small  values  of  Ra,  Ma  or  Pe  indicate  that 
the  thermal  convection  within  the  melt  is  of  a  minor  effect.  Hence,  neither 
thermocapillary  nor  buoyancy  induced  convection  is  included  in  the  model;  forced 
convection  is  included  here  mainly  to  account  for  the  mass  conservation 
constraint.   Accordingly,  in  the  present  model,  besides  the  transient  and  two- 
dimensional  conductional  effects  within  the  melt  and  solid,  the  only  convection 
effect  considered  is  the  forced  flow  fluid  caused  by  the  continuous  pulling  of  the 
crystal  at  the  speed  Up. 

4.5  Mathematical  Model 
In  the  following,  a  theoretical  model  of  heat  transfer  in  the  EFG  process 
under  transient  conditions  is  formulated.   The  governing  equations  of  the  heat 
transfer  process  in  the  melt  and  solid  are  presented  along  with  the  constraint 
controlling  the  shape  of  the  melt/solid  interface.  The  interface  condition  adopted 
here  is  based  on  phase  change  equilibrium  (isothermal  phase  change  condition) 
and  local  heat  balance  without  considering  the  microscopic  effects  such  as  the 
Gibbs-Thompson  effect  or  the  anisotropic  surface  energy  (Pelce  1988). 
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We  are  interested  in  the  axisymmetric  situation,  i.e.  r-z  coordinates,  with 
which  the  crystal  is  grown  from  the  melt.  A  schematic  of  the  melt  region  of  the 
present  EFG  model  along  with  the  associated  boundary  conditions  is  given  in 
Figure  4-4.  Table  4-2  gives  the  physical  dimensions  and  the  material  properties  of 
sapphire,  while  Table  4-3  summarizes  the  important  parameters  and  their 
numerical  values  encountered  in  the  present  model  based  on  Table  4-2.  Both 
large  St  and  small  St  cases  are  numerically  investigated.  The  values  of  Bond 
number,  Rayleigh  number,  Marangoni  number,  and  Peclet  number  are  all  small 
here,  and  hence  conduction  is  the  dominant  heat  transfer  mechanism  within  the 
melt  and  the  crystal,  and  capillarity  dominates  static  pressure  in  terms  of 
controlling  the  shape  of  the  meniscus.  Under  such  a  circumstance,  the  meniscus 
profile  is  close  to  a  straight  line  (Kalejs  et  al.  1983a).  Regarding  the  governing 
equations  and  boundary  conditions  relevant  to  the  present  system,  the  following 
presents  the  Laplace- Young  equation  (controlling  the  shape  of  the  melt/gas 
interface),  energy  transport  equation,  mass  continuity  equation,  and  movement  of 
the  solid/melt  interface  in  the  r-z  axisymmetric  coordinate  system. 


(i)  meniscus  shape 

The  Laplace- Young  equation  governing  the  equilibrium  meniscus  shape  in  the 
axisymmetric  geometry  is 
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where  7,^  is  the  surface  tension  between  the  melt  and  the  air,  f  describes  the 
meniscus  shape,  and  Ap  is  the  difference  in  applied  pressure  between  the  liquid 
phase  and  the  gas  phase.  The  bottom  edge  of  the  meniscus  is  defined  by  the  die 
radius,  r^,  ,  and  the  top  edge  is  defined  by  the  crystal  radius,  r,  . 
(ii)  energy  equation 
Melt 

The  energy  equation  in  the  liquid  phase  in  dimensional  form  is 

dt     -  bz        r  dr      dr      dz  dz 
where  the  subscript  1  denotes  the  liquidus  phase  condition,     is  the  convection 
speed  along  the  axial  direction  and  f  and  h  define  the  meniscus  shape  in  r  and  z 
direction,  respectively.  The  boundary  conditions  are  given  as  follows  : 

(a)  central  line 

.^=0,       r=0,  Oizi/i(0)  ("^-22) 
dr 

(b)  the  melt  inlet 

r,  =  r,.       0^r^r,,z=0  (4-23) 

where  Tb  is  the  base  temperature  at  the  melt  inlet,  and  rb  is  the  radius  of  the 
contact  line  between  the  meniscus  and  die,  and 

(c)  along  the  side  surface  of  the  melt,  heat  is  transferred  to  the  surroundings  by 
the  combined  conduction,  convection  and  radiation 
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dT 

-k'  =  hfT^  -r )  +  eo[{T^  r  -r  %       r=M,  0^z<h^  (4-24) 
an 

where  n  is  normal  unit  vector  on  the  meniscus  surface,  h,  the  heat  transfer 
coefficient,  e  the  emissivity,  a  Stefan-Boltzman  constant,  T,  the  surroundings 
temperature,  and  h,  the  height  of  trij  unction  point. 
Crystal 

Similarly,  the  energy  equation  in  the  solid  phase  is 

pCp(^^u—l)=±l.(kr-^)+±(k—i),     Q^r:Lr,h<z<l  (4-25) 
'    '  dt     '  dt        rdr  '  dr      Bz  '  dz 

where  the  subscript  s  denotes  the  solidus  phase  condition,  r,  is  the  radius  of  the 

trijunction  point,  and  1  is  the  height  of  the  domain  modelled.   As  to  the  boundary 

conditions, 

(a)  along  the  side  surface  of  the  solid,  heat  is  transferred  to  the  surroundings  by 
the  combined  conduction,  convection  and  radiation 

dT 

-kr^  =  KiT-L)  ^  ^<^iiTr-V],       r-r{z)MO^z<l  (4-26) 

(b)  along  the  top  boundary  of  the  solid,  the  heat  flux  is  specified 

PsCPs^T,  =      ,       0<r<r(/),z  =  /  (4-27) 

dz 

where  G,  is  the  prescribed  heat  flux  out  of  the  top  boundary  from  the  crystal,  and 

(c)  along  the  centerline  a  symmetric  condition  is  applied 
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dT 

 i=0,       r=0,  /j(0)izi/  ("^-28) 

dr 

(iii)  mass  conservation 

Because  of  the  low  Peclet  number,  convection  does  not  affect  heat  transfer 
in  a  significant  way.  The  velocity  in  the  melt  is  determined  using  the  one- 
dimensional  continuity  equation: 

where  v^is  the  advection  speed  in  the  melt,  which  is  determined  by  the  local  radius. 

(iv)  conditions  at  solid/melt  interface 

Finally,  at  the  melt/sold  interface  z  =  h(r,t),  both  equilibrium  solidification 
and  conservation  of  energy  are  satisfied 

T  =  J  =  T  (4.30  a) 

Ism 


-  kf-Il  -  pMu- ^\  (4-30  b) 

'  dn     '  dn       '  ^   '    dt  " 


The  various  issues  of  the  specification  of  an  appropriate  contact  condition  at 
the  trijunction  point  have  been  discussed  in  Chapter  3.  Conventionally,  a  constant 
contact  angle,  in  the  spirit  of  the  Young's  equilibrium  contact  condition,  has  been 
used  (Sackinger  et  al.  1989,  Derby  and  Brown  1986,  Ettouney  and  Brown  1983, 
Boley  1975),  namely, 


-S  =  (u  -'-^)tan(0(r)  -  <^^  (4-31) 
where  r,  is  the  radius  of  the  trijunction  point,  h,  is  the  height  of  the  trijunction 
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point,  (^(t)  is  the  instantaneous  contact  angle  of  the  meniscus  at  the  trijunction 
point,  (i>o  is  the  equilibrium  contact  angle. 

The  above  equations  are  then  written  in  dimensionless  form,  by  using  the 
following  reference  length,  time,  and  temperature,  respectively 
[Length]:       1^  =  h  (radius  of  the  die), 
[Time]:         t,  =  h^/a^,  (diffusion  time  scale), 
[Temperature]: AT  =  T^  -  T„ 

Using  the  upper  case  letters  to  designate  the  nondimensional  parameters  in 
space,  T*  =  t  /t^  to  designate  nondimensional  time,  and 
T;  =     /  a,  ,  i  =  i,  s 
p  =  p,  /p, ,     h  =  h,  /  h„   k  =  k,  /  k,, 
Up  =\i,tJ  a,  V,  =  V,  rb  /  a, 

=  (T,  -  TJ  /  AT,  9,  =  (T,  -  TJ  /  AT,  G,*  =  G,  /  (p,Cp,ATa,  /  rj 

The  dimensionless  energy  equation  and  corresponding  boundary  conditions 
for  the  liquid  phase  become 
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=  l±iR^)^±i^),     0<R^FiZ),  0<Z<H{R)  (4.32) 
dr'       dZ       R  dR     dR      dZ  dZ 


—L=0,  R=0,  OiZi//(0)  ("^-33) 
dR 

0=1,         0  ii?  i/?  ,Z  =0  (4-34) 

-  ^  =  Bi/e,  -  ej  -  Rad^lid,-^*  -  (6,^-^1       J?  =  F(Z),  0  ^  Z  ^  i/,  (4.35) 

where  R,  Z,  N,  F,  and  H  are  all  nondimensionalized  with  respect  to  r^  . 

The  dimensionless  energy  equation  and  the  corresponding  boundary 
conditions  for  the  solid  phase  become 


36       36  \  3      96       3  36 

(^+U^)  =  T[1—{R^)  +_l(_i)],    0<R<R{Z),  0<Z<L  (4-36) 
dr'     'dZ  R  3R     dR      3Z  3Z 


^  =5jA-e.)  -M[(e,*^'  -(Q.^^V    R=RiZ),H^^-Z^L  (4.37) 


36 

U6  -r_L  =  G  ,       Oii?i/?(/)  ,  Z  =L  (4.38) 
"  '  3Z 


36 

—1=0,  R=0,  //(O)iZiL  (4-39) 
3R 

where  Bii  =  h,  rb  /  k,,  Bi,  =  h,  rb  /  k,  are  the  Riot  number,  and 

Rad,  =  ea  rb  AT^  /  k,  ,  Rad,  =  ea  r^,  AT^  /  k,  are  the  radiation  number. 
At  the  interface  S(R,Z,r)  =  Z  -  H(R,r)  =  0,  the  temperature  equilibrium 

and  local  heat  balance  become 
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e  =  0  (4.40  a) 


where  St  is  the  Stefan  number. 

The  interface  shape  is  determined  by  thermal  equilibrium  and  local  heat 
balance.  The  exact  location  of  the  trij  unction  point  which  defines  the  radius  of 
grown  crystal  is  decided  by  enforcing  a  constant  contact  angle  at  the  upper 
trijunction  point,  which  has  been  employed  by  several  authors  (Sackinger  et  al. 
1989,  Derby  and  Brown  1986,  Ettouney  and  Brown  1983,  Boley  1975),  i.e., 

^  =  (t;  -—1)  tan(0(r)  -<^^  (4.41) 

where  Up  is  the  pulling  speed,  <^(t)  is  the  instantaneous  contact  angle  of  the 
meniscus  at  the  trijunction  point,  and  <^o  the  contact  angle  of  the  meniscus  at 
steady-state.   At  the  lower  trijunction  point,  on  the  other  hand,  the  contact  radius 
is  fixed  by  the  sharp  edge  of  the  die  tip. 

Rewriting  the  heat  equation  in  general  form  for  both  liquid  and  solid 
phases,  ^i(r,z,t),  i  =  l,s: 

(— +v;.— )  =  r.[-L-r_(/?_i)  +_Z_(^)]  (4.42) 
ar*    *az      '  /?a/?   m    dz  dz 

where  V,  =  V,  and  V,  =  Up  .   Employing  the  transformation  from  the  (R,Z,r*) 
coordinates  to  curvilinear  the  (^,17,7)  coordinates, 
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R  =  R(^rj,T)  , 


e  =  |(R,Z,r*), 


Z  =  Z(?,ij,t)  , 


7j  =  rj  (R,Z,r'), 


(4.43) 


T 


T 


=  T 


the  above  general  equation  after  dropping  the  subscript  "i"  can  be  rewritten  as 


where  q,  =  Z,^  +  R,^       q^  =  Z^Z,  +  R^R„  q,  =        +  Rj', 

J  =  Z^R,  -  Z,R^ ,      W„  =  V,R,  ,  and      W,  =  -  V.R^ 
Regrouping  the  above  equation,  and  treating  the  grid  movement  induced 
convection  terms  and  cross  derivative  terms  as  pseudo-source  terms,  we  move 
them  to  the  right  hand  side  of  the  equation: 


where  the  spatial  derivative  terms  are  treated  based  on  the  second-order  central 
difference  scheme  in  the  context  of  second-order  finite  volume  formulation,  and 
the  unsteady  term  is  by  the  backward  Euler  method. 

As  already  discussed,  the  computational  strategy  designed  here  utilizes  an 
explicit  front  tracking  technique,  based  on  a  combined  Lagrangian  and  Eulerian 
approach,   to  locate  and  advance  the  moving  boundaries  between  the  crystal  and 
the  melt,  as  well  as  the  new  meniscus  shape.  In  this  procedure,  the  markers  along 


(4.44) 


(4.45) 
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the  interface  are  used  to  define  and  advance  the  location  of  the  interface  in  time; 
the  movement  is  tracked  via  a  Lagrangian  framework.  The  velocity  of  each 
marker  is  determined  by  the  balance  between  the  difference  of  thermal  gradients 
in  the  solid  and  in  the  melt,  and  the  latent  heat  absorption  as  given  in  Eq. 
(4.40  b).  Based  on  the  newly  defined  boundary  shapes,  a  grid  will  be  generated 
again  and  the  rate  of  grid  movement  evaluated  accordingly.  The  energy  equation 
in  each  phase  will  be  computed  in  an  Eulerian  framework.  All  these  procedures 
are  performed  in  a  fully  coupled  manner  involving  interactions  among  the 
temperature  field,  meniscus  and  interface  motion,  and  grid  movement  at  each 
iteration.   Within  each  time  step,  the  iterations  continue  until  all  the  governing 
equations,  boundary  conditions,  and  interface  constraints  are  satisfied  within  a 
specified  convergence  criterion;  here,  it  requires  that  the  summation  over  the 
domain  of  the  absolute  imbalance  of  the  nondimensional  fluxes  in  each 
computational  cell  to  be  smaller  5xlO"^ 


4.6  Numerical  Simulation  of  the  EFG  Process 
The  grid  distribution,  generated  by  an  algebraic  procedure,  transfinite 
interpolation  technique  (Thompson  et  al.  1985),  is  clustered  in  the  regions  close 
to  the  central  symmetric  line  and  the  outside  boundary  of  the  crystal  and  melt. 
The  control  surfaces  are  determined  with  a  power-law  clustering  scheme  over  half 
of  the  computational  domain 
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=  ^^(4r^'"^     i  =0,l,...,iVr  (4.46  a) 

2  N 

r 

2.  =  !!^{Lr,  j=o,u...,Nz  (4.46b) 

where  2Nr  +  1  and  2N^  +  1  are  the  total  number  of  control  volumes  along  R- 
and  Z-  direction,  respectively.  The  exponents  m  and  n  control  the  rate  of 
clustering  in  the  R-  and  Z-  direction,  respectively.  A  value  of  1.55  for  m  and  1 
for  n  is  used  in  this  study,  i.e.  the  clustering  is  assigned  along  the  R-  direction 
only.  The  grid  distribution  is  symmetric  about  the  lines  R  =  F(R,t)/2  and  Z  = 
H(R,r)/2.   The  grid  with  finer  resolution  in  the  region  close  to  the  outside 
boundary  of  the  crystal  and  melt  is  found  important  for  yielding  good  solution 
accuracy. 

The  complete  determination  of  the  meniscus  profile  involves  the  solution 
of  the  Laplace- Young  equation  constrained  by  hydrostatic  and  externally 
enforced  pressure  difference,  Ap,  the  Bond  number,  defined  in  Eq.  (4.20),  the 
meniscus  contact  angle  at  the  trij  unction  point,       the  meniscus  height  of  the 
trijunction  point  H„  and  the  die  radius,  R^.    There  are  also  issues  of  static  versus 
dynamic  contact  angles,  as  examined  in  Chapter  3.  One  set  of  solufions  obtained 
using  the  linear  meniscus  profile,  i.e., in  non-dimensional  terms,  R^  =  1,  H,  = 
0.188,R,  =  0.763,  is  chosen  to  compare  the  solution  of  Laplace- Young  equation 
of  Bo  from  2.8x10-^0  2.8xlO-\with  Ap/Apgr,  =  5%  and  0%.  The  computed 
meniscus  profiles  are  quite  close  to  linear,  as  shown  in  Figure  4-5.  With  the 
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straight  line  assumption,  the  meniscus  volume  does  not  differ  by  more  than  2  % 
from  that  calculated  based  on  the  Laplace- Young  equation.    This  observation  is 
consistent  with  the  assumption  made  in  Kalejs  et  al.  (1983a).  For  these 
parameters  relevant  to  the  sapphire  fibre  production,  and  under  the  condition  that 
the  static  contact  angle  controls  the  trijunction  point  movement,  the  use  of  a 
straight  meniscus  profile  is  a  reasonable  approximation.    Both  large  St  (Case  1:  St 
=  1)  and  small  St  (Case  2:  St  =  0.024)  are  computed.   Thermophysical  properties 
and  physical  dimensions  of  sapphire  (AI2O3),  processing  parameters  and  the 
corresponding  non-dimensional  operating  parameters  are  given  in  Tables  4-2  and 
4-3.  As  can  be  seen,  for  simplicity,  constant  properties  are  used  for  both  solid 
and  liquid  phases.  This  aspect  is  not  expected  to  affect  the  solution  in  a  major 
way. 

Due  to  the  symmetry  of  the  problem,  only  half  of  the  physical  domain  is 
modelled.   For  the  entire  physical  domain  encompassing  both  the  crystal  and  the 
melt,  effectively  a  total  of  41x41  nonuniform  grid  is  employed.   Due  to  the 
complex  nature  of  the  problem,  such  a  grid  size  already  requires  substantial 
computing  resources  to  support  the  numerical  simulations.   In  the  present  work, 
both  engineering  workstations,  including  IRIS  Indigo,  DEC  5000  and  3100,  and 
the  Cray  Y-MP  supercomputer  have  been  used.  For  those  workstations,  it 
typically  requires  ten  hours  or  more  of  CPU  time  to  complete  a  simulation. 

The  initial  condition  of  both  St  =  1  and  St  =  0.024  is  shown  in  Figure  4-6, 
and  the  corresponding  steady-state  solution  is  shown  in  Figure  4-7  and  4-8, 
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respectively.    Density  and  thermal  diffusivity  are  taken  to  be  the  same  for  both 
the  melt  and  solid  phases.  Figure  4-7  shows  the  isotherms  of  the  steady-state 
solutions  of  St  =  1  in  both  phases,  along  with  the  grid  lines  used  to  conduct  the 
computation.    The  interface  is  flat  in  the  inner  region,  and  becomes  more  curved 
toward  the  outer  region  due  to  the  heat  loss  to  the  ambient  environment.    It  is 
noted  that  because  the  outer  boundary  of  the  solid  and  melt  are  of  different 
slopes,  the  local  curvatures  of  the  isotherms  in  these  two  phases  are  different. 

The  grid  system  and  the  steady-state  isotherm  distribution  of  St  =  0.024 
are  shown  in  Figure  4-8;  due  to  the  small  value  of  Stefan  number,  the  thermal 
gradients  in  the  melt  and  the  crystal  regions  are  of  different  magnitudes. 
Furthermore,  it  is  interesting  to  note  that  with  St  =  0.024  the  solid/melt  interface 
is  convex  toward  the  melt,  while  in  previous  section,  it  is  found  that  with  St  =  1 
the  interface  is  convex  toward  the  solid  region. 

In  terms  of  the  computational  practices,  we  discuss  two  scaling  procedures 
utilized  to  non-dimensionalize  and,  more  critically,  normalize  the  various  terms 
present  in  the  governing  equations  and  their  implications  on  the  computational 
efficiency.  As  it  turns  out,  care  should  be  taken  in  choosing  the  appropriate 
reference  quantities  because  for  many  phase  change  problems  this  choice  has  a 
major  impact  on  the  amount  of  computing  cost  required  to  conduct  a  simulation. 
It  is  noted  that  in  the  present  problem,  the  two  energy  transport  equations  are 
applied  in  the  regions  separated  by  a  moving  interface,  causing  a  nonlinear 
coupling  between  the  two  equations.   Hence,  in  the  course  of  solving  these  energy 
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equations  with  the  backward  Euler  time  stepping  scheme,  while  a  von  Neumann 
type  of  stability  analysis  indicates  that  the  computation  is  unconditionally  stable 
for  each  of  the  energy  equations,  the  nonlinear  coupling  caused  by  the  interface 
movement  makes  it  only  conditionally  stable.  Numerical  experiments  indicate 
that,  as  expected,  the  range  of  the  time  step  size  acceptable  for  a  stable 
computation  depends  on  the  distribution  of  the  grid  points  in  each  phase. 

Table  4-4  presents  two  different  scaling  procedures  and  the  resulting  non- 
dimensional  equations  of  energy  transport  and  interface  movement.   As  can  be 
seen,  the  only  difference  between  the  two  procedures  is  the  velocity  scale;  in 
choice  1,  a  standard  characteristic  diffiisional  velocity  is  used,  while  in  choice  2, 
the  Stefan  number  is  included  in  addition  to  the  diffusional  velocity  scale.  The 
main  motivation  of  the  second  choice  stems  from  the  observation  that  as  can  be 
easily  deduced  from  Eq.  (4.40  b),  with  choice  1,  the  non-dimensional  speed  of 
interface  movement  is  scaled  with  St.  Accordingly,  in  non-dimensional  terms,  for 
the  low  St  cases,  the  interface  moves  at  a  correspondingly  slow  rate.  As  already 
mentioned,  even  with  the  use  of  an  implicit  procedure,  the  computation  is  not 
unconditionally  stable.  Numerical  experiments  have  indicated  that,  in  terms  of 
the  non-dimensional  values,  the  time  step  sizes  of  the  set  of  governing  equations 
resulting  from  both  choice  1  and  2  have  very  comparable  stability  restrictions. 
Hence,  considering  the  fact  that  the  time  scale  of  choice  2  is  larger  than  that  of 
choice  1  by  a  factor  of  1/St,  the  relative  computing  efficiency  of  the  two  scaling 
procedures  depends  highly  on  the  value  of  the  Stefan  number.   As  an  illustration. 
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computations  have  been  conducted  using  both  scaling  procedures  and  starting  with 
the  identical  initial  conditions  for  the  low  St  case  studied  here.  In  both 
computations,  the  non-dimensional  time  step  size  is  At  =  10'^  ,  which  is  about  10 
times  of  the  nondimensional  diffusional  time  scale  based  on  the  smallest  grid 
spacing,  namely  (Ax)^Va  ,  where  (Ax)^  is  the  smallest  grid  spacing  and  a  is  the 
thermal  diffusivity.  For  both  scaling  procedures,  numerical  experiments  indicate 
that  this  value  of  At  is  very  close  to  the  stability  limit  of  the  implicit  Euler  time 
stepping  scheme  adopted  in  the  present  work. 

Figure  4-9  compares  the  paths  of  the  computing  histories  with  both  scaling 
procedures  from  an  identical  initial  condition  toward  the  steady-state  solution 
subject  to  the  same  boundary  conditions.  The  Stefan  number  chosen  in  the 
present  work  is  0.024.  Substantially  different  convergent  behaviors  have  been 
observed  between  the  two  scaling  procedures.   As  already  mentioned,  the 
differences  of  the  convergence  characteristic  is  largely  caused  by  the  different 
orders  of  magnitude  of  the  non-dimensional  interface  velocity  yielded  by  the  two 
scaling  procedures.   With  choice  1,  the  interface  speed  is  of  the  order  of  St,  while 
the  choice  2,  it  is  of  the  order  of  1.  Consequently,  choice  2  needs  a  far  fewer 
number  of  non-dimensional  time  steps  to  reach  the  steady  state  solution.  With 
choice  1,  however,  the  system  appears  to  go  through  a  long  transient  period 
before  it  can  reach  the  steady-state.   As  a  check,  the  converged  steady-state 
solution  obtained  from  choice  2  is  substituted  into  the  computation  based  on 
choice  l,and  this  solution  is  immediately  accepted  as  the  converged  solution, 
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confirming  that  with  both  scaling  procedures,  the  computations  eventually  lead  to 
the  same  steady-state  solution.  However,  the  computing  costs  of  the  two  scaling 
procedures  differs  by  an  order  of  (1/St).   Obviously,  choice  2  is  appropriate  for 
the  low  St  computations. 

4.7  Conclusions 

In  this  Chapter,  a  thermocapillary  model  has  been  developed  to  simulate 
the  characteristics  of  the  EFG  process.  The  model  is  axisym metric,  solving  the 
energy  equation  in  both  solid  and  melt  phases  separately,  and  tracks  the 
movement  and  evolutions  of  their  interface  explicitly  using  a  combined 
Lagrangian/Eulerian    method.  The  meniscus  behavior  is  governed  by  the  Laplace- 
Young  equation  subject  to  a  specified  contact  angle  at  the  trijunction  point. 

It  has  been  established  that  for  the  thin  fibre  process,  the  value  of  various 
control  parameters  such  as  the  Stefan  number,  Rayleigh  number,  Marangoni 
number,  and  Bond  number  are  modest,  and  hence  conduction  is  the  primary  heat 
transfer  mechanism  within  the  melt  and  the  crystal,  and  capillarity  dominates 
static  pressure  in  terms  of  controlling  the  shape  of  the  meniscus. 

The  quasi-stationary  approximation  for  interface  is  found  to  be  inaccurate 
for  capturing  the  phase  boundary  except  for  small  Stefan  numbers.   With  a  full 
treatment,  the  temperature  field,  interface  calculation  and  grid  movement  all 
interact  over  each  interaction  and  time  step.  Such  a  computational  procedure  is 
applicable  to  an  isothermal  interface.  A  generalized  interface  motion  technique  is 
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developed,  which  is  not  limited  to  a  single-value  or  isothermal  interface. 
Computations  show  that  this  method  is  both  accurate  and  versatile. 

Comparison  of  the  two  scaling  procedures  for  the  interface  movement 
shows  that  the  characteristic  velocity  scale  should  be  defined  to  ensure  that  the 
solid/liquid  interface  advances  at  a  nondimensional  speed  of  order  1 .  Since  the 
interface  movement  is  primarily  responsible  for  the  nonlinearity  of  the  phase 
change  problem,  a  conventional  choice  based  on  the  diffusion  scale,  choice  1,  will 
yield  a  set  of  equations  that  requires  the  computing  cost  to  increase  by  a  factor  of 
1/St. 

It  is  also  interesting  to  note  that  with  St  =  1 ,  the  solid/liquid  interface  is 
convex  toward  the  solid  region,  while  with  St  =  0.024,  the  interface  is  convex 
toward  the  melt. 
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Table  4-1  Assessment  of  time  stepping  schemes  for  phase  change  problem. 

(a)     Comparisons  of  the  predicted  interface  location  at  r  =  0.2  using  backward 
Euler  and  Crank-Nicolson  time  stepping  schemes  for  Eqs.  (4.3)  and  (4.4) 
with  11  and  21  grids. 


11  grid           (Ax  =  0.0158) 

21  grid           (Ax  =  0.0079) 

At 

Backward  Euler 

Crank-Nicolson 

Backward  Euler 

Crank-Nicolson 

10-^ 

0.223 

0.222 

0.227 

0.223 

10-^ 

0.226 

0.225 

0.226 

0.226 

10^ 

0.226 

0.226 

0.226 

0.226 

10-^ 

0.226 

0.226 

0.226 

0.226 

(b)     Comparisons  of  the  required  CPU  time  on  Micro Vax  3100  using 

backward  Euler  and  Crank-Nicolson  time  stepping  schemes  for  Eqs.  (4.3) 
and  (4.4)  with  11  and  21  grids. 


11  grid           (Ax  =0.0158) 

21  grid           (Ax  =  0.0079) 

At 

Backward  Euler 

Crank-Nicolson 

Backward  Euler 

Crank-Nicolson 

10-^ 

0:12.28 

0:13.07 

0:16.77 

0:21.00 

10-^ 

0:39.84 

1:12.24 

1:09.79 

1:13.68 

10^ 

4:54.87 

4:5967 

9:15.15 

9:14.34 

10-^ 

49:49.27 

49:48.20 

1:31:27.93 

1:31:15.27 
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Table  4-2      Thermophysical  properties  and  physical  dimensions  of  sapphire 
(AI2O3)  used  in  the  present  EFG  simulation  (Case  1:  St  =  1  and 
Case  2:  St  =  0.024). 


The  material  property  of  sapphire  (ALO^')  Case  1  Case  2 

(St  =  1)  (St  =  0.024) 

k,  (W/cm-K)  0.1  0.1  0.1 

p,  (g/cm')  3.05  3.05  3.05 

Cp,  (J/g-K)  1.26  1.26  1.26 

e,  0.9  0.9  0.9 

k,  (W/cm-K)  0.1  0.1  0.1 

p,  (gW)  4.00  3.05  3.05 

Cp,  (J/g-K)  1.26  1.26  1.26 

e,  0.9  0.9 

T„  (K)  2316  2316  2316 

lf(J/g)  1046  25  1046 

<t>o  (deg)*  135  135 

h(W/cm2-K)  5  1.1x10-^ 

h  (cm)  0.02  0.02 

AT(K)  20  20 

Ta  (K)  2216  316 

Up  (cm/s)  0.13  0.062 

(1/K)  3x10-^  3x10-^ 

7  (dyne/cm)  0.69  0.69 

a  (W/cm    -  K")  5.67x10  '^      5.67x10  '^ 

H  (g/cm-s)  0.59  0.59 

f(cm-/s)  0.19  0.19 


'  (t>„  is  defined  as  the  angle  between  solid/melt  interface  and  melt/gas  interface. 
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Table  4-3      Definition  and  magnitude  of  key  dimensioniess  parameters  based  on 
values  given  in  Table  4-1  (Case  1:  St  =  1  and  Case  2:  St  =  0.024). 


Case  1:  St  =  1 

Case  2:  St  =  0.024 

a.t 

Fo  =  —  :       Fourier  number 

2 

0(1) 

0(1) 

St  =  —       ;       Stefan  number 

I 

J 

1 

0.024 

Bo  =              ;       Bond  number 

2.8x10^ 

2.8x10-* 

Ra  =                ;      Rayleigh  number 

v^ 

4.2x10"^ 

4.2x10'^ 

1        1  Arr, 
Ma  =  ;  Marangoni  number 

1.56 

1.56 

Pe  =         ;       Peclet  number 
«/ 

0.1 

0.05 

hr. 

Bi  =  — —  ;       Biot  number 
k 

1 

with  0^  =  -5 

2.2x10-^ 
with  0^  =  -100 

Rod  =   ■  ;    Radiation  number  . 

K 

0 

1.8x10' 
==^=-.  ^  

135 

Table  4-4      Two  different  scaling  procedure  and  resulting  nondimensional 
energy  equations 

(a)  choice  1 

length  scale:  1,  =  r^  ,         velocity  scale:     =  a,  /  1^  ,  time  scale:  t^  =  Ir  /     ,  and 

temperature    scale:  AT  =  difference  between  temperature  at  melting  point  and  at  die  tip. 

ae.     dd.    1  a   ae,    a  ae, 

energy  eq,  for  melt    —  +  F  —  =  -—(R—i)  +  — (— ) 

dx      'dZ      RdR    dR      dZ  dZ 


ae,     ae,       13    ae,    a  ae, 

ensery      for  solid    ^  *  t/,^  -  [||(«^)  ^ 


^/  D  a/f 

interface  movement  eg.  — -  -k — -   =   -t- (1/  - — )„ 

^        dN       dN        St  ^    dx  ^ 

where  R  =  r/1,  ,  Z  =  z/1,  ,  Up  =  Up/v,  ,  V,  =  v,/v,  ,  r  =  t/  t,  ,  0,  =  (T.-TJ/AT, 

=  (T,-TJ/AT,  p  =  p,  /  p,  ,  and  a  =  a,  /  a, 

(b)  choice  2 

length  scale:  1,  =  rb  ,  velocity  scale:  v,  =  St  a,  /  1,  , 

time  scale:  t,  =  1^  /  v^ ,  and  temperature  scale:  AT 

.  at    ^az     /^a/?  br    bz  dz 


energy  eq.  for  solid    St(^  +  F-^)  =  a  [1-^{R^)  +  — (-^)] 

ax      'dZ  RdR     dR      dZ  dZ 


interface  movement  eq.  — -  -  k — -   =   o  (U  -  — 

aiv     aiv         p  dx  ^ 
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Table  4-4  -  continued 
Boundary  conditions  based  on  choice  1 

(1)  for  the  melt: 

dQ, 

symmetric  b.c.  :  — -=Q      @R  =  Q,  OiZ iH(Q) 

dR 

inlet  temperature  :       6,  =  1 ,       @  0  i    i/?^ ,  Z  =  0 


80  jT  T  ' 

convectivelrcu&ative  b.c.  :     "  ^  =  Bif,e,  -     +  ^^(6,+-^*  -  (e^*-^*] ,     @  R  =  FiZ) ,  0  i  Z  nH^ 


(2)  for  the  solid: 


dd 

symmetric  b.c.  :       —  =0  ,       @  R  =  0,  //(O)  <Z  <L 

dR 


89 

convecti\el diffusive  b.c:    U^B^-V—  =  G^,     ®Q<R:iR{L)  ,  Z=L 


50  T  T 

convectivelradiative  b.c.  :    -  — 7  =  flj,(e,  -  6.)  +  RadjiiQ^*-^*  -  (O,*-^*] ,    ®R=R{Z),  HiZsL 


(3)  at  the  interface  S(R,Z,7)  =  Z  -  H(R,r)  =  0  : 

isothermal  b.c.  :  9=0 


local  heat  balance  :     [(V5)  J'  =  A  (f/  -  HL) 

'      St     "     dr  " 


0 
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meit  1 

solid 

\^  j 

%  =0 

S(X,t) 


Governing  equations: 
energy  eq.  for  melt 


dd  _  d^Q  ^ 

5t      dX^  dY^ 


_  50  ^  p^dS 
dN  ~  St  dx 


interface  movement  eq.       - —  =  -i— 


B.C. 

8  =  1,  X=Q,  T 

6=0,  XnS,  t 


I.e. 


e(x,t  =0.1)  =  1 


er/[X) 

5(X.T,=0.1)  =  2Ay^ 


Neumann's  solution: 


6  =  1- 


erf^X) 


S  =  2X^ 


ke^'erf{,X)  =  ^ 


Figure  4-1  Exact  solution  for  melting  of  a  solid  confined 
region. 


to  a  semi-infinite 
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Figure  4-2  Comparison  of  interface  trajectories  predicted  by  the  coupled 
(full  equations)  and  decoupled  (quasi-stationary  equations) 
approaches  with  three  values  of  the  Stefan  numbers. 


139 


Figure  4-3  Comparison  of  computed  and  exact  solutions  for  St  =  0.1303 
St  =  1.2216  and  St  =  2.8576.  The  figures  on  the  left  are  the  ' 
exact  and  computed  temperature  fields  at  different  time 
instants.   The  figures  on  the  right  show  the  superposed  exact 
and  computed  interface  locations  versus  time. 


■1 
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(a)  schematic  configuration  in  the  melt  region 


Qs  Cp,  U,  T,  -  k,^  =  Gs 


(b)  boundary  condtions 


Figure  4-4  Schematic  configurations  and  boundary  conditions  of  tiie 
present  EFG  process. 
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(a)  Ap  =  5  % 


1^ 


0.73  0.8  0.83  0.9  0.93  1 


f/rb 

(b)  Ap  =  0  % 


Figure  4-5  Comparison  of  linear  meniscus  profile  with  the  Laplace- Young 

solution  with  R,  =  I,  H,  =  0.188.  =  0.763,  and  various  values 
of  Bond  number,    (a)  Ap  =  5  %  and  (b)  Ap  =  0  %. 
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(a)  grid  system 


Figure  4-6  J^;^^  ^^i^.,.r:n.^  co„.ou.  of  U,e  inicia,  co„di.io„s 
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(a)  grid  system 


(b)  isotherm  (A0  =  0.2) 


Figure  4-7  The  grid  system  and  isothermal  contours  of  the  steadv-state 
solution  for  Case  1:  St  =  1.  ' 
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(a)  grid  system 


(b)  isotherm  (A0  =  0.2) 


Figure  4-8  The  grid  system  and  isothermal  contours  of  the  steady-state 
solution  for  Case  2;  St  =  0.024. 
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CHAPTER  5 

DYNAMIC  SIMULATION  OF  THE  EFG  PROCESS 


5.1  Introduction 

In  the  EFG  technique,  as  in  other  meniscus-defined  growth  methods,  the 
processes  of  heat  and  mass  transfer  serve  as  the  link  between  the  crystal  growth 
condition  and  its  resulting  shape  and  structure.   For  example,  the  radius  of  the 
crystal  is  controlled  by  the  interaction  of  heat  transfer  from  the  melt  to  the  crystal 
and  the  ambient,  with  the  shape  of  melt/gas  meniscus  controlled  by  the  force 
balance  between  surface  tension,  the  pressure  difference  across  the  interface,  and 
the  local  shear  stress  distribution.   Accordingly,  the  radius  of  the  crystal  can  be 
affected  by  either  the  rate  of  heat  transfer  from  the  system  to  the  surroundings,  by 
varying  the  temperature  of  the  melt  that  enters  the  die,  or  by  changing  the  crystal 
pulling  rate. 

In  practice,  both  the  die  tip  temperature  and  the  puller  speed  may 
experience  either  intentional  or  unintentional  variations  in  time.  The 
unintentional  variations  can  result  from  the  perturbations  to  the  operating 
conditions  caused  by  environmental  fluctuations,  coil  power  fluctuations,  or  motor 
speed  irregularities.   The  intentional  variations  can  be  designed  to  compensate  the 
aforementioned  perturbations  to  achieve  a  controlled  growth  process  and  uniform 
crystal  properties  and  dimensions.   In  order  to  minimize  the  variations  in  fibre 
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quality  caused  by  these  fluctuations  in  manufacturing  conditions  some  type  of 
active  control  scheme  may  be  necessary  (Gevelber  et  al.  1988,  1987,  Gevelber  and 
Stephanopoulos   1987,  Meric  1979  and  Robinson  1971).  Furthermore,  so  far  in 
the  literature,  the  static  contact  condition  is  normally  adopted.   Effort  has  been 
made  to  develop  two  different  models  for  the  contact  condition  at  the 
solid/melt/gas   trijunction  point;  they  are:  (1)  the  conventional  static  model  which 
requires  that  the  Young's  equilibrium  contact  condition,  namely,  a  fixed  contact 
angle  between  the  surfaces  separating  melt/gas  and  solid/melt  be  maintained 
(Crowley  1983,  Ettouney  and  Brown  1983,  and  Dussan  V.  1979),  and  (2)  a 
dynamic  model  which,  instead  of  fixing  the  contact  angle,  allows  it  to  vary 
according  to  the  magnitude  and  direction  of  the  instantaneous  movement  of  the 
trijunction  point. 

The  purpose  of  this  chapter  is  to  mathematically  analyze  the  transport 
phenomena  in  the  EFG  system  so  as  to  establish  the  connections  between  the 
actual  growth  conditions  and  the  corresponding  solidification  characteristics.  Two 
different  values  of  Stefan  number,  St  =  1  and  0.024,  have  been  computed  for  the 
simulation.   For  St  =  1,  both  temperature  fluctuations  at  the  die  tip  and  speed 
fluctuations  of  the  puller  motor  are  considered;  for  St  =  0.024,  the  pull  speed 
fluctuations  are  investigated.   Also  addressed  in  the  St  =  0.024  case  is  the 
behavior  of  the  trijunction  point  constrained  by  a  dynamic  contact  condition. 
Relevant  issues  identified  and  insight  gained  from  the  simulation  are  discussed. 
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5.2  Dynamic  Perturbation  on  the  EFG  Process  with  St  =  1 
The  steady-state  solution  for  St  =  1  shown  in  Figure  4-7  is  used  as  the 
initial  condition  for  the  dynamic  simulations.  Two  types  of  external  perturbations 
will  be  considered:  the  fluctuation  of  temperature  and,  accordingly,  the  incoming 
heat  flux,  at  the  lower  boundary  of  the  melt,  and  the  fluctuation  of  pull  speed. 
The  fluctuations  will  be  represented  either  by  combinations  of  sinusoidal  waves  or 
by  a  saw-toothed  profile. 

5.2.1  Base  Temperature  Perturbation 

We  first  consider  the  dynamic  characteristics  of  the  present  EFG  process  in 
response  to  a  base  temperature  perturbation.    Varying  the  inlet  temperature 
changes  the  amount  of  heat  flux  into  the  melt  which,  in  turn,  affects  the 
solidification  rate,  the  meniscus  height  H(R),  and  through  the  contact  angle 
constraint,  the  trij  unction  location  R,.  In  this  section  we  examine  the  effect  of  the 
base  temperature  perturbation  by  keeping  the  pull  speed  (Up  =  0.1,Pe  =0.1) 
and  all  other  operating  parameters  unchanged.    The  disturbance  of  the  base 
temperature  is  given  as 

^(^)  =(U.,-i;^..sin(^)  (5.1) 
.-1  Wf 

where  (0b)avg  is  the  average  base  temperature  ,  i.e.  (6^),^^  =  1, 

Aj  is  the  amplitude  of  the  disturbance  of  the  i*  harmonic,  and 
Ui  is  the  period  of  the  disturbance  of  the  i""  harmonic. 
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A  time  step  size  of  At  =  10'^  is  used  for  all  the  calculations.   Two  forms 
of  perturbation  will  be  considered  here;  one  with  a  single  harmonic  and  another 
with  three  harmonics.  For  the  perturbation  with  a  single  harmonic,  the  period,  fi, 
is  either  500  At  or  20  At,  both  with  the  same  amplitude,  A  =0.1.  For  the 
perturbation  with  three  harmonics,  the  periods  of  the  three  harmonics  are, 
respectively  500  At  ,  292  At  and  108  At,  each  with  the  same  magnitude  of  0.1/3. 
In  both  cases,  as  shown  in  Figures  5-1  and  5-2,  clear  correlations  exist  between  . 
He  and  R,.  Firstly,  while  H,  closely  follows  the  movement  of  d^,  there  is  a  phase 
angle  between  them.  The  phase  angle  is  20  At  (14.4°)  for  a  single  harmonic 
perturbation  of  the  period  of  500  At,  8  At  (144°)  for  a  single  harmonic  of  the 
period  of  20  At,  and  12  At  (8.6°)  for  the  case  of  three-harmonic  perturbation. 
The  phase  angle  between      and  6^,  indicates  the  time  required  for  the 
temperature  perturbation  at  the  base  of  the  meniscus  to  travel  to  the  melt/solid 
interface. 

Figure  5-1  also  reveals  that  for  both  cases,      and  Rj  are  out  of  phase. 
This  feature  is  a  direct  outcome  of  the  constraint  imposed  at  the  trij unction  point. 
As  indicated  by  Eq.  (4.41),  with  a  constant  Up  and  a  fixed  contact  angle  between 
the  solid/gas  and  melt/gas  interfaces,  the  rate  of  change  of      and  R^  are  of  the 
opposite  signs.  Another  aspect  worth  mentioning  is  the  noniinearity  of  the  system. 
Although  both  H,  and  R,  move  at  the  identical  frequencies  as  those  of  their 
time-averaged  mean  values  are  different  from  the  steady-state  values.  This 
phenomenon  indicates  that  there  are  accumulated  effects  caused  by  the  nonlinear 
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physics  intrinsic  to  the  process.  Finally,  Figure  5-1  also  clearly  demonstrates  that 
there  is  a  large  difference  in  the  system's  sensitivity,  in  terms  of  the  magnitudes  of 
He  and  Re  variations,  to  low  and  high  frequency  forcings. 

The  interface  shapes  and  locations  of  the  melt/solid  at  selected  time 
instants  are  plotted  in  Figure  5-3  for  both  the  single-harmonic  and  the  three- 
harmonic  forcings.  While  the  maximum  magnitudes  of  the  temperature 
fluctuations  in  both  cases  are  the  same,  the  ranges  of  the  solid/melt  interface 
movement  are  not  identical.  The  case  with  a  single-harmonic  perturbation 
(Q=500At)  exhibits  more  interface  movement  than  that  with  three-harmonic 
perturbations.    This  phenomenon  largely  reflects  the  fact  that  each  harmonic 
affects  the  interface  movement  in  a  different  time  scale.  It  should  be  noted  that 
due  to  the  multi-dimensional  nature  of  the  problem,  the  heat  flux  difference 
across  the  solid/melt  interface  is  not  uniform,  resulting  in  varying  interface 
speeds,  as  evidenced  by  the  nonuniform  distances  between  the  interfaces  shown  at 
different  time  instants  shown  in  Figure  5-3. 

To  better  understand  the  effect  of  the  forcing  frequency  of  6^,  on  the 
solid/melt  interface  movement,  the  phase  angles  among  He  ,  Up  ,  and  are 
inspected.   By  keeping  the  amplitude  of  the  size  of  fluctuation  unchanged,  i.e. 
A =0.1,  and  varying  the  period  from  500  At  to  20  At,  the  computed  results  are 
summarized  in  Figure  5-4.  Figure  5-4  (a)  shows  that  the  phase  shift  between  He 
and  Up  increases  as  the  frequency  of  the  disturbance  increases.  He  and  Re  are 
always  out  of  phase  because  the  rates  of  change  of  He  and  Re  are  of  the  opposite 


151 

signs  when  Up  is  fixed,  as  the  trijunction  constraint,  Eq,  (4.41),  implies.  The  small 
phase  angle  variation  between      and      is  caused  by  the  solid/ melt  interface 
curvature  at  the  trijunction  point.  The  percentages  variation  of  H,.  and  with 
respect  to  the  forcing  frequency  are  presented  in  Figure  5-4  (b).  Both  the 
variations  of  H,  and      decrease  as  the  forcing  frequency  increases.  However, 
the  characteristics  of     and  R^  variations  behave  differently  with  respect  to  the 
forcing  frequencies.   For  example,  for  the  time  period  of  &b  perturbation  larger 
than  200  At,  the  solidification  process  responds  to  the  external  forcing  with 
approximately  the  same  sensitivity.  As  the  forcing  frequency  becomes  sufficiently 
high,  however,  the  sensitivities  of  both  H,  and  R,  decrease  noticeably,  indicating 
that  the  thermal  inertia  causes  the  system  not  to  be  able  to  keep  pace  with  the 
high  forcing  frequencies. 

5.2.2 Pull  Speed  Perturbation 

The  pull  speed  in  the  typical  EFG  process  (2-8  cm/min)  is  one  to  two 
orders  of  magnitude  larger  than  that  of  conventional  growth  methods,  e.g. 
Czochraiski  (Zulehner  1983)  or  floating-zone  (Duranceau  and  Brown  1986,  and 
Ostrach  1980,  1983).  For  the  process  considered  here,  the  radius  of  a  grown 
crystal  typically  is  only  a  few  microns.  As  already  discussed,  convection  is  not 
important  for  the  thermal  transport  in  either  the  melt  or  the  crystal.  However, 
variation  in  pull  speed  has  a  direct  impact  on  the  solid/melt  interface  and  the 
trijunction  point  movements,  as  reflected  by  Eq.  (4.41).  Increasing  pull  speed 
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affects  both  the  amount  of  heat  transported  by  advection  and  the  amount  of 
latent  heat  released  upon  solidification,  which  cause  the  meniscus  height  to 
increase  and  the  solid/melt  interface  to  deform. 

For  the  results  to  be  presented  in  the  following,  \]^{t)  is  perturbed  while 
all  other  parameters  are  unchanged.   The  steady-state  condition  represented  in 
Figure  4-7  is  used  as  the  initial  condition  to  start  the  computations.    Two  types  of 
perturbation  are  considered,  namely,  the  sinusoidal  and  the  saw-toothed  shapes. 
With  the  sinusoidal  perturbations,  Up(T)  varies  as 


where  (Up),vg  =  0. 1  is  average  pull  speed.  This  choice  of  the  shape  of 
perturbations  is  motivated  by  the  desire  to  investigate  the  system  dynamics  with 
respect  to  (1)  few  discrete  harmonics,  and  (2)  a  collection  of  harmonics  with  a 
wide  distribution  of  frequencies. 

We  first  study  two  cases  of  a  single-harmonic  perturbation  with  J]  =  500  At 
and  20  At,  respectively,  both  with  the  same  fluctuating  magnitude,  A  =  0.5.  The 
time  histories  of  Up  ,  H,  and      are  shown  in  Figure  5-5.  While  both  H,  and  R, 
respond  to  the  Up  perturbation  at  the  same  frequency,  which  is  similar  to  the  6^, 
forcing,  the  decrease  in  the  system's  sensitivity  with  respect  to  forcing  frequency 


(5.2) 


and  with  the  saw-toothed  perturbation,  Up  (t)  varies  as 


U,ir)  =  (UX,      rim)  forO<r<  fi/2, 

(^p)ayg  ■''^ (r-Q/2)  /(Q /2)  -A    forW2  <T<Q 


(5.3) 
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seems  more  pronounced  for  U,  than  for  6,  .  There  is  also  a  clearly  identifiable 
phase  angle  between  H,  and      ,  which  for  the  period  of  500  Ar  is  45  Ar  (32.4^) 
and  for  20  Ar  is  5  Ar  (90°).  This  phase  angle  represents  the  time  scale  needed  by 
the  thermal  field  in  the  melt  to  respond  to  the  perturbation  imposed  by       ;  it  is 
larger  than  that  caused  by  the  base  temperature  fluctuations  with  a  low  frequency 
perturbation,  but  smaller  with  a  high  frequency  perturbation.    The  other 
noticeable  aspect  is  that  H,  and  R,  are  no  longer  out  of  phase,  which  is  different 
from  the  case  of  d,  perturbations.    Again,  the  reason  for  this  phenomenon  can  be 
found  from  the  constraint  at  the  trijunction  point.  Previously,  Up  was  a  constant, 
resulting  in  the  requirement  that  dH,/dr  and  dR,/dr  are  essentially  of  opposite 
signs.  Here,  Up  is  no  longer  unchanged,  and  accordingly,  the  relationship  between 
dH,/dr  and  dR,/dr  is  no  longer  unchanged.   The  phase  angle  between  the  time 
histories  of  H,  and      is  a  reflection  of  this  characteristic.    By  keeping  the  forcing 
amplitude  unchanged,  i.e., A  =  0.5,  and  varying  the  period  from  500  Ar  to  20  Ar, 
the  relative  phase  angles  between  H,  and  Up,  between  H,  and  K  ,  as  well  as  the 
percentage  variations  of  H,  and  R„  are  depicted  in  Figure  5-6.  Figure  5-6(a) 
shows  that  the  phase  shift  between  H,  and  Up  increases  as  the  forcing  frequency 
increases.  The  phase  shifts  between  both  H,  and  Up  ,  and  between  H,  and  R, 
become  larger  as  the  frequency  becomes  higher.  Figure  5-6(a)  indicates  that  H, 
and  R,  tend  to  be  more  out  of  phase  as  the  forcing  frequency  increases.  The 
trend  of  the  phase  angle  between  H,  and  external  forcing  is  largely  the  same,  as 
evidenced  by  the  similarity  between  -Figures  5-4(a)  and  5-6(a).  The  variations  of 
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H,  and      with  respect  to  the  forcing  frequency  are  presented  in  Figure  5-6(b). 
Sensitivities  of  both  H,  and  R,  decrease  as  the  forcing  frequency  increases.  It 
indicates  that  by  perturbing  either  6^,  or  Up,  the  thermal  inertia  of  the  system 
cannot  respond  equally  with  respect  to  different  harmonics.   However,  unlike  the 
case  of  By,  forcing  shown  in  Figure  5-3(b),  Figure  5-6(b)  indicates  that  the 
sensitivities  of  both  H,  and  R,  reduce  more  or  less  at  the  same  rates  as  the 
forcing  frequency  of  Up  increases. 

We  next  study  the  movement  of  the  solid/melt  interface  in  response  to  the 
perturbed  Up  containing  three  incommensurate  harmonics,  with  the  periods  of  500 
At,  292  At  and  108  At,  all  with  the  same  magnitude  of  0.5/3.  The  frequencies 
imposed  here  are  identical  to  a  case  of  perturbing  6^,  studied  previously,  shown  in 
Figure  5-2.  Figure  5-7  shows  the  time  histories  of  Up  ,  H,  and  R,  for  this  case. 
While  Up  in  Figure  5-7  and  6^  in  Figure  5-2  are  of  the  identical  form,  the  shapes 
and  magnitudes  of  variations  in  H,  and  R,  are  quite  different  in  the  two  cases. 
With  the  0b  perturbation,  all  three  forcing  frequencies  are  better  captured  by  both 
He  and  R^  .  As  already  observed,  by  perturbing  Up  ,  both  H,  and  R,  progressively 
reduce  their  sensitivities  with  respect  to  the  external  forcing;  by  penurbing  ^^„  on 
the  other  hand,  H,,  and  especially  R„  maintain  their  sensitivities  with  respect  to 
the  forcing  frequency  until  it  becomes  quite  high. 

To  further  illustrate  this  characteristic,  a  saw-toothed  perturbation  is  given 
to  both  Up  and  6^,  with  a  period  of  20  At  and  A  =0.1.  The  time  histories  of  Up, 
He  and  R^  are  shown  in  Figure  5-8.  With  such  high  frequency  perturbations,  while 
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the  system  responds  to  the  forcing  in  a  modest  manner,  the  sensitivity  to  Up  is 
much  smaller  than  to  • 

5.3  Dynamic  Perturbation  on  the  EFG  Process  with  St  =  0.024 
Next,  we  choose  a  Stefan  number  of  0.024,  which  is  more  representative  of 
the  current  state-of-the-art  in  actual  production.  Together  with  previous  section 
St  =  l,the  dynamic  characteristics  of  the  EFG  process  under  different  conditions 
can  be  assessed  and  compared.   For  the  results  presented  here,  all  the  relevant 
conditions  are  given  in  Table  4-3.  Compared  to  the  case  of  St  =  1,  the  steady- 
state  pull  speed  is  reduced  by  half  and  some  modification  are  made  to  the 
thermal  boundary  conditions.  The  motivation  for  making  these  changes  is  to 
bring  the  modeled  condition  closer  to  a  typical  production  condition.  The  most 
important  change  is  the  Stefan  number,  other  modification  do  not  change  the 
qualitative  features  presented  below.  The  simulation  starts  with  the  steady-state 
solution.  Figure  5-8,  as  the  initial  condition.  Figure  5-9  shows  the  time  series 
plots  of  two  cases,  one  with  a  fluctuating  pull  speed  of  the  period  of  500At, 
At=10''  ,and  another  with  20  At.  Shown  there  are  the  time  histories  of  the  pull 
speed.  Up  ,  the  height  and  the  radius  of  the  trijunction  point,  H,  and  R,  , 
respectively.  Consistent  with  the  case  of  St  =  1  studied  in  previous  section,  the 
crystal  diameter  responds  to  the  externally  imposed  perturbation  at  the 
corresponding  frequency,  but  the  sensitivity  of  the  response  depends  on  the 
frequency  of  perturbation.    For  a  slower  perturbation  of  the  period  of  500At,  both 
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H,  and  R,  exhibit  higher  levels  of  oscillation  than  for  a  faster  perturbation  of 
20At;  the  differences  between  them  are  substantial.   Furthermore,  also  consistent 
with  the  case  of  St  =  I,  the  nonlinearity  of  the  underlining  physics  has  caused  the 
time  averaged  values  of  H,  and  R,  to  be  different  from  the  steady-stale  values. 

Figure  5-10  gives  an  overall  depiction  of  the  percentage  variations  of  H, 
and  R,  with  respect  to  different  perturbation  periods  predicted  by  the  present 
model.  For  both      and  R,  ,  the  tendency  of  exhibiting  reduced  sensitivities  to 
the  external  perturbations  as  the  frequency  increases  can  be  clearly  observed. 
The  reason  for  this  phenomenon  obviously  lies  in  the  relative  competition 
between  the  time  scale  of  the  grower  and  the  time  scale  of  the  perturbation. 
Since  At  is  10'^  ,  for  the  pull  speed  fluctuation  of  the  period  of  several  hundred 
At,  the  system's  dominant  non-dimensional  time  scale,  which  is  of  order  1 
according  to  choice  2  given  in  Table  4-3,  and  the  perturbations  time  scale  are 
comparable.   Hence,  the  system  is  able  to  respond  and  follow  the  perturbation 
closely,  resulting  in  a  clear  connection  between  the  dynamic  behaviors  of  Up  and 
H,  and  R^  .  For  high  frequency  perturbations,  of  the  period  of  hundred  At  or 
less,  the  system's  time  scale  is  too  long,  causing  decreases  in  sensitivity  of  and 
Rj  ,  as  observed  in  Figure  5-10.  This  phenomenon  has  also  been  observed  in  the 
actual  growth  system.  Figure  5-11  shows  a  Bode  diagram  based  on  a  series  of 
experiments  to  depict  the  frequency  response  between  the  fibre  diameter  and  pull 
speed  (Backman  1992).  It  indicates  that  with  relatively  slow  varying  perturbations 
over  a  range  of  frequencies  in  pull  speed,  the  fibre  diameter  responds  with  a 
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certain  sensitivity;  as  the  perturbation  frequency  increases  beyond  a  threshold,  the 
system  can  no  longer  respond  to  it  with  equal  sensitivities.  Hence,  the  predicted 
trend  of  the  system  response  is  qualitatively  verified  experimentally. 

5.4  Dynamic  Contact  Angle 
While  it  seems  reasonable  to  use  the  Laplace- Young  equation  to  determine 
the  meniscus  shape  due  to  the  smallness  of  Ra,  Pe  and  Ma,  it  is  an  appropriate 
prescription  of  the  contact  condition  at  trijunction  point.  Two  models  have  been 
investigated  in  this  regard,  including  a  conventionally  adopted  static  model  which 
fixes  the  contact  angle,  and  a  dynamic  model  which  regulates  the  contact  angle 
according  to  the  speed  and  direction  of  the  instantaneous  movement  of  the 
trijunction  point.  The  primary  focus  of  this  section  is  to  investigate  the  effect  of 
the  perturbation  in  pull  speed,  subject  to  two  different  models  for  the  contact 
condition  at  trijunction  point,  on  the  characteristics  of  the  dimensional  variation 
of  the  fibre. 

We  present  the  results  obtained  based  on  a  dynamic  model  for  the  contact 
condition  at  the  trijunction  point.  The  model,  as  schematically  depicted  in  Figure 
5-12,  now  allows  different  values  of  contact  angle  according  to  the  instantaneous 
direction  as  well  as  the  speed  of  movement  of  the  trijunction  point.  Drawing  the 
analogy  to  the  advancing  and  receding  contact  angle  of  a  liquid  drop  on  a  flat 
surface  (Dussan  V.  1979),  two  different  contact  angles  are  assigned  as  the 
asymptotic  values  for  the  outward  and  inward  moving  cases,  respectively. 
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Furthermore,  if  the  contact  point  moves  at  a  slow  enough  speed,  in  the  range  of 
+e  and  -e  as  shown  in  Figure  5-12,  then  the  contact  angle  is  allowed  to  vary 
between  the  two  asymptotic  angles.  For  the  case  implemented,  we  take  the 
asymptotic  value  of  the  contact  angle  associated  with  the  outward  moving  case  to 
be  the  same  as  that  adopted  in  the  static  model,  namely,  135°.  The  asymptotic 
contact  angle  of  the  inward  moving  case  is  taken  as  134°,  and  e  is  taken  as  0.1.  In 
other  words,  the  variation  of  contact  angle  between  the  inward  and  outward 
moving  cases  is  1°.  While  this  variation  seems  very  small,  it  has  a  quite  noticeable 
impact  on  the  characteristics  of  the  solidification  process. 

Figure  5-13  shows  the  comparison  of  the  time  histories  of  both  H,  and 
between  the  cases  with  static  and  dynamic  contact  conditions,  where  a 
perturbation  of  Up  with  a  period  of  100  At  is  enforced.  It  is  striking  to  observe 
that  with  a  seemingly  small  modification  of  the  contact  condition,  the  resulting 
values  of  H,  and  R,  are  substantially  affected.  For  R,  ,  both  time-averaged  and 
fluctuating  values  are  quite  different  with  the  two  contact  conditions;  the  dynamic 
contact  model  yields  a  crystal  of  a  somewhat  larger  mean  diameter  and  smaller 
variations.   For      ,  on  the  other  hand,  the  time-averaged  fibre  diameter  is 
smaller  with  the  use  of  the  dynamic  contact  condition;  the  fluctuating  magnitudes 
in  both  cases,  however,  appear  quite  comparable.    Hence,  there  seems  to  be 
differences  in  behavior  between  H,  and  R,  resulting  from  the  static  and  dynamic 
models.   In  this  regard,  several  points  can  be  noted  to  elucidate  the  results 
observed  here.  As  observed  in  Figure  4-8,  the  solid/melt  interface  is  concave 
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toward  the  solid  phase,  indicating  that  as  the  contact  angle  is  reduced,  the 
trijunction  point  tends  to  move  outward,  yielding  a  larger  diameter  as  shown  in 
Figure  5-13.  This  will  also  seem  to  cause  a  corresponding  increase  in  H,  ,  as  seen 
in  the  initial  transient  behavior  in  the  time  series  plot  of      .  However,  a  larger 
crystal  diameter  also  causes  a  larger  volume  and  a  higher  heat  loss  both  via 
conduction  and  advection.   Accordingly,  there  is  also  a  change  in  the  balance  of 
the  thermal  constraint  at  the  interface,  which  eventually  reduces  the  value  of  H, 
as  observed  in  the  later  stage  of  Figure  5-13. 

It  should  be  noted  that  if  the  solid/melt  interface  were  of  the  opposite 
curvature,  then  the  trend  of  R,  will  change  as  well.  Hence,  the  interaction  of 
thermal  and  contact  condition  plays  a  pivotal  role  in  determining  the  solidification 
characteristics.    Furthermore,  since  the  contact  angle  is  no  longer  fixed  in  the 
dynamic  model,  an  extra  degree  of  nonlinearity  also  appears  via  the  enforcement 
of  the  contact  condition.  Plots  of  the  rate  of  change  of  H,  and  R,  are  given  in 
Figure  5-14  to  demonstrate  this  point.  Clearly,  with  the  static  contact  condition,  a 
single  harmonic  reflecting  that  of  the  imposed  perturbation  of  Up  is  reproduced  in 
both  H,  and  R,  .  On  the  other  hand,  with  the  dynamic  contact  condition,  the 
movement  of  trijunction  point  is  no  longer  consistent  with  the  imposed  frequency, 
instead,  it  exhibits  a  quasi-periodic  pattern  with  many  additional  frequencies 
present.  These  additional  frequencies  appear  because  the  way  the  contact  angle, 
and  hence  the  values  of  H,  and  R,  ,  vary  in  response  to  the  pull  speed 
perturbations. 
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In  general,  the  value  of  the  equilibrium  contact  angle  may  be  difficult  to 
quantify.  Here  an  assessment  is  made  to  investigate  the  effect  of  contact  angle  on 
the  solidification  dynamics.  Three  values  of  static  contact  angle,  namely,  107°, 
125°  and  135°,  are  used.  Figure  5-15  shows  the  frequency  responses  of  solid/melt 
interface  with  different  static  contact  angles,  to  the  pull  speed  fluctuations.  The 
overall  features  of  the  interface  movement  of  the  three  contact  angles  show 
similar  trend  for  both  AH,  and  AR,  ,  although  the  variations  of  AH,  and  AR,  are 
affected  by  the  value  of  the  contact  angle.  Because  the  curvature  of  the  interface 
near  the  trij  unction  point  is  sensitive  to  the  local  energy  balance  and  does  not 
vary  monotonically  with  the  contact  angle,  there  is  no  simple  trend  between  AH,  , 
ARj  and  <i)^  . 

5.5  Summarv  and  Concluding  Remarks 
In  this  Chapter,  a  thermocapillary  dynamic  model  has  been  utilized  to 
simulate  the  EFG  process.  The  dynamic  behavior  of  the  fibre  growth  process 
with  respect  to  the  external  disturbance,  including  the  base  temperature 
perturbation  and  the  pull  speed  perturbation,  is  investigated.   Based  on  the  results 
presented,  the  following  conclusions  are  made: 

l.For  both  base  temperature  and  pull  speed  perturbations,  the  system 
reaches  the  limit  cycle  with  the  same  frequency  as  the  external  forcing  after  a 
period  of  transient  response;  however,  the  magnitudes  of  their  variations  with 
respect  to  these  perturbations  are  found  to  be  dependent  on  the  forcing 
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frequencies.    It  has  also  been  predicted  that  as  the  forcing  frequency  increases, 
the  solid/melt  interface  responds  at  a  reduced  sensitivity.  This  trend  has  been 
observed  experimentally.   The  time-averaged  values  of      and      are  different 
from  the  steady-state  values  because  of  the  nonlinearity  of  the  fibre  growth 
process. 

2.  While  the  Stefan  number  substantially  affects  the  temperature 
distribution,  interface  locations,  and  crystal  diameter,  it  does  not  affect  the 
qualitative  behavior  of  the  solidification  process  in  the  dynamic  environment. 

3.  The  phase  shifts  as  well  as  magnitudes  of  H,.  and      are  found  to  be 
forcing  frequency  dependent  for  both  base  temperature  and  pull  speed 
perturbations,    (i)  The  phase  angles  between  H,  and  external  perturbation 
increase  as  the  forcing  frequency  increases,  (ii)  Both  H,  and  R,  are  less  affected 
by  the  high  frequency  harmonics  of  the  perturbations.    However,  with  the  Up 
perturbation,  responses  in  H,.  and  R^  decrease  at  relatively  uniform  rates  as  the 
forcing  frequency  increases,  but  with  the     perturbation,        and  R,  exhibit  two 
different  sensitivity  scales  with  respect  to  the  forcing  frequency. 

4.  The  relationship  between  H,  and  R,  is  dictated  by  the  constraint  at  the 
trijunction  point,  which  differs  in  two  different  types  of  forcing  mechanisms 
studied.  With  the  base  temperature  perturbation,  while      and  R,  remains 
essentially  out  of  phase  at  all  frequencies.  H,  and     maintain  a  phase  angle  which 
is  frequency  dependent.    With  the  pull  speed  penurbation,  the  relationship 
between  He  and  R,  ,  as  well  as  between  H^  and  Up  ,  is  frequency  dependent. 


162 

5.  Two  different  models  have  been  utilized  for  the  contact  condition  at 
trij unction  point,  including  a  conventionally  adopted  static  contact  angle,  and  a 
dynamic  model  which  allows  the  contact  angle  to  vary  according  to  the  direction 
and  the  rate  of  the  trij  unction  point  movement.   It  turns  out  that  the  two  models 
can  yield  substantial  differences  in  the  crystal  diameter,  including  the  time- 
averaged  as  well  as  the  fluctuating  magnitudes.   Furthermore,  additional 
harmonics  can  also  be  yielded  by  the  dynamic  contact  condition  due  to  the  extra 
nonlinearity  contained  by  it.  These  observations  have  clearly  indicated  that  the 
static  contact  condition  is  not  sufficient  to  yield  all  the  information  regarding  the 
dynamic  behavior  of  the  crystal  correctly. 

6.  Available  experimental  information  has  confirmed  the  characteristics  of 
frequency  response  of     discussed  here. 
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Figure  5-1  Time  histories  of  H,  and  R,  subject  to  a  single  harmonic 
perturbation  of     with  (a)  Q  =  500  Ar  and  (b)  Q  =  20  At 
(St  =  1.0) 


Figure  5-2  Time  histories  of      and  R,  subject  to  a  three-harmonic 
perturbation  of     with  12,  =  500  Ar,  ^  =  292  At,  and 
=  108  Ar  .  (St  =  1.0) 
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(b)  a  three-harmonic  pemirbauon  in  6^  with  Q,  =  500  At,      ~  292  At  and  Q, 
=  108  At  . 

Figure  5-3  The  interface  shapes  at  different  time  instants  of  (a)  a  single- 
harmonic  perturbation  in  d^,  with  Q,  =  500  At-  and  (b)  a  three- 
harmonic  perturbation  in  6^,  with  Q,  =  500  Ar,      =  292  Ar,  and 
=  108  Ar  .  (St  =  1.0) 
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Figure  5-5  Time  histories  H.  and       subject  to  a  single  harmonic 
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ure  5-7  Time  histories  H,  and      subject  to  a  three-harmonic 

perturbation  of  with  fi,  =  500  A.,  =  292  Ar  and 
"3  -  108  At  .  (St  =  1.0) 
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Figure  5-9  Time  histories  of      and      subject  to  a  single  harmonic 
perturbation  of  U,  with  (a)  period=500  Ar  and 
(b)  period =20  Ar.  (St  =  0.024) 
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Figure  5-13  Time  histories  of  H,  and      subject  to  a  single  harmonic 
perturbation  of  Up  with  period  =  lOOAr  .  (St  =  0.024) 
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Figure  5-15    Characteristics  of  interface  movement  with  respect  to  different  forcing 
periods  of  Up  for  different  static  contact  angles  (St  =  0.024) 


CHAPTER  6 
CONCLUSIONS  AND  RECOMMENDATIONS 


This  chapter  presents  some  conclusions  that  have  been  drawn  from  this 
work,  together  with  a  few  suggestions  as  to  how  this  work  could  be  usefully 
extended. 

6.1  Achievements  and  Findings 

An  adequate  theoretical  formulation  of  the  EFG  process  needs  to  address 
several  issues.  They  include  the  method  of  predicting  the  thermal  fields  in  the 
solid  and  the  liquid,  the  shapes  of  solid/liquid  and  liquid/gas  interfaces,  front 
tracking  and  moving  grid  techniques,  understanding  of  the  meniscus  behaviors  and 
the  contact  condition  at  the  trij unction  point,  and  appropriate  scaling  procedures. 
In  this  study,  a  thermocapillary  model  has  been  developed  to  investigate  the 
dynamic  characteristics  and  associated  transport  phenomena  of  the  EFG  process, 
subject  to  the  external  perturbations  and  different  contact  conditions.   Based  on 
the  computed  results,  the  following  conclusions  are  made. 

An  order  of  magnitude  estimation  shows  that  the  convective  mechanisms 
are  minor  influences  on  the  transport  process  in  a  typical  EFG  process; 
conduction  dominates  energy  transport.   In  the  absence  of  an  externally  imposed 
pressure  field,  capillarity  dominates  the  shape  of  the  meniscus. 
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The  shape  of  the  liquid/gas  meniscus  that  connects  the  growing  crystal  with 
the  die  tip,  governed  by  the  Laplace- Young  equation  and  subject  to  the  contact 
condition  at  the  trijunction  point,  is  critical  in  determining  the  properties  of  the 
solidified  crystal.  A  study  of  meniscus  formation  shows  that  multiple  solutions 
may  exist.  The  realizability  of  these  possible  meniscus  profiles  has  been  examined 
using  the  minimum  free  energy  concept.  Interactions  among  Bond  number, 
pressurization,  aspect  ratio,  contact  condition,  and  pulling  direction  have  been 
found  to  result  in  a  variety  of  meniscus  behaviors.  A  hysteresis  model  for  the 
dynamic  contact  line  motion  has  been  employed  to  investigate  the  quasi- 
equilibrium  dynamics  of  the  meniscus.  It  is  found  that  the  relation  in  terms  of 
waveform  and  phase  of  the  radius  of  the  trijunction  point  with  respect  to  the 
height  of  the  trijunction  point  cannot  be  fully  captured  by  the  formula 
conventionally  employed  in  crystal  growth  simulations.   Full  fluid-dynamic 
consideration  is  required  to  assess  the  implication  of  such  behavior. 

A  new  generalized  interface  tracking  technique  is  developed,  which 
overcomes  restrictions  of  single-valuedness  of  the  interface  imposed  by  commonly 
used  mapping  methods.  This  method  is  shown  to  be  both  accurate  and  versatile. 
Combining  this  method  with  a  moving  grid  technique,  detailed  thermal  fields  in 
both  solid  and  liquid  phases  can  be  computed. 

Two  different  scaling  procedures  for  the  interface  movement  are  employed 
for  both  St  =  1  and  St  =  0.024.  One  scaling  procedure  uses  the  diffusional 
characteristic  velocity  as  the  velocity  scale,  and  the  other  includes  both  the 
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diffusional  velocity  and  the  latent  heat  effects.  Since  the  interface  movement  is 
primarily  responsible  for  the  nonlinearity  of  the  phase  change  problem,  the 
procedure  which  accounts  for  the  latent  heat  effect  results  in  better  computational 
efficiency  for  low  St  cases. 

The  sensitivities  and  dynamic  characteristics  of  the  EFG  system  subject  to 
the  base  temperature  and  pull  speed  perturbations  are  explored.   For  both 
disturbances,  the  crystal  radius      and  the  height  of  the  trijunction  point  H,  reach 
limit  cycles  at  the  same  frequency  as  the  external  forcing;  the  magnitude  of  their 
variations  with  respect  to  these  perturbations  are  found  to  be  dependent  on  the 
forcing  frequencies.   As  the  forcing  frequency  increases,  the  solid/liquid  interface 
responds  with  a  reduced  sensitivity.  This  trend  has  been  observed  experimentally. 
Also,  the  time-averaged  values  of     and  R,  are  different  from  the  steady-state 
values  because  of  the  nonlinearity  of  the  solidification  process. 

Two  different  models  have  been  utilized  to  simulate  the  contact  condition 
at  the  trijunction  point,  including  a  conventionally  adopted  static  contact 
condition,  and  a  dynamic  model  in  which  the  contact  angle  varies  according  to  the 
direction  and  the  simultaneous  speed  of  the  trijunction  movement.   It  turns  out 
that  the  two  contact  conditions  can  yield  substantial  differences  in  solidification 
characteristics  as  reflected  by  the  variation  in  crystal  diameter.   These  results 
indicate  that  the  conventional  static  model  is  not  adequate  to  model  some  aspects 
of  the  trijunction  point  behavior,  particularly  in  the  presence  of  the  phase  change. 
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6.2  Directions  for  Future  Research 

The  present  thermocapillary  model  has  produced  some  interesting  results, 
and  has  been  partially  confirmed  experimentally.   It  will  offer  useful  insights  into 
a  practically  important  process.  To  further  extend  the  content  of  the  model, 
further  work  can  be  conducted  in  research  areas  as  suggested  below. 

The  contact  condition  specified  at  the  trij  unction  point  is  found  to  be 
important  in  determining  the  crystal  growth  characteristics.    Our  simplified 
dynamic  contact  model  is  a  first  step  toward  accounting  for  the  nonequilibrium 
effects  at  the  trij  unction  point.  Hydrodynamic  problems  involving  dynamic  contact 
lines  are  unsteady,  highly  nonlinear,  and  involve  free  and  moving  interfaces  with 
singular  solution  near  the  moving  contact  lines  (Shopov  and  Bazhlekov  1991, 
Shopov  et  al.  1990,  Dussan  V.  1979).  An  adequate  theory  pertaining  to  such 
systems  is  still  in  the  formative  stages,  and  much  work  remains  to  be  done. 
Experimental  information  is  critical  in  helping  to  gain  a  better  understanding  of 
this  issue. 

In  this  study,  the  phase  boundary  is  determined  by  the  condition  of  local 
thermodynamic  equilibrium.   This  assumption  is  only  valid  for  small  undercooling 
and  slow  pull  speed.  In  general,  the  interface  temperature  is  a  function  of  the 
local  rate  of  phase  change  kinetics,  introducing  an  extra  degree  of  freedom  which 
is  constrained  by  a  kinetic  relationship  between  interface  velocity  and 
undercooling  (Wang  and  Matthys  1992,  Clyne  1984,  Levi  and  Mehrabian  1982, 
Shingu  and  Ozaki,  1975).  A  complete  model  capable  of  predicting  accurately  the 
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solid/liquid  interface  velocity  and  interface  location  according  to  both  the  thermal 
and  kinetic  constraints  will  be  useful. 

On  a  practical  side,  the  findings  reported  here  can  be  used  to  facilitate 
design  of  appropriate  active  control  schemes  to  improve  the  fibre  quality 
(Gevelber  et  al.  1988,  1987,  Gevelber  and  Stephanopoulos   1987,  Meric  1979, 
Robison  1971). 
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